Rによる二次分析入門

「子どもの生活と学びに関する親子調査」を用いて

はじめに

0.1 本コースの概要

本コースでは，Rをこれまで使ったことのない方，あるいは他のソフトウェアからRへの移行を検討している方を対象に，RStudioを用いた二次分析の基礎を習得することを目指す．

具体的には，以下の内容を扱う．

• プロジェクト管理の方法
• SSJDAから提供されるデータの読み込みと加工
• 記述統計量の算出と可視化
• 回帰分析の実行と解釈
• パネルデータの構築と基礎的な分析手法

0.2 使用するデータ

データは，ベネッセ教育総合研究所よりSSJデータアーカイブに寄託されている「子どもの生活と学びに関する親子調査」関連の3つを使用する．

• メインデータ：Wave1からWave7までの2015年から2021年にかけて収集されたデータ（調査番号：1571）
• 補助データ1：2016年度実施の語彙力・読解力調査（調査番号：1577）
• 補助データ2：2019年度実施の同調査（調査番号：1578）

0.2.1 データの概要

項目	内容
調査期間	2015年〜2021年（7時点）
対象	小学1年生〜高校3年生の子どもと保護者
サンプル規模	約20,000組（有効回答は各Wave約15,000〜16,000組）
データ形式	.sav, .csv, .dta

0.2.2 変数の構成

変数は以下の順で並んでいる．

1. フェイス変数（学年，性別，父母学歴，父母職業）
2. Wave1〜Wave7の子ども調査データ
3. Wave1〜Wave7の保護者調査データ

0.2.3 変数名の命名規則

データセット内の変数名は，以下のルールに従って命名されている．

パネル調査項目（複数時点で同じ質問をたずねている変数）

• 形式：dsXXXXXXXXXX_wYc または dsXXXXXXXXXX_wYp
• ds + 10桁の数字 = 変数ID（同じ質問は全Waveで同じIDを持つ）
• _w + 数字 = Wave番号（1〜7）
• c = 子ども調査（child） / p = 保護者調査（parent）

例：

変数名	意味
ds0000000080_w1c	Wave1・子ども調査「朝起きる時間」
ds0000000080_w2c	Wave2・子ども調査「朝起きる時間」（同じ質問）
ds0000000080_w1p	Wave1・保護者調査の同様の質問

変数IDの活用

同じ質問は全Waveで同じ変数ID（10桁の数字部分）を持つ．複数時点のデータを扱う際に，この規則を利用して変数を特定できる．

派生変数・作成変数

データ提供元で作成された派生変数は，異なる命名規則を使用している．

変数名	意味
PFED15〜PFED21	父親の最終学歴（15=2015年, 21=2021年）
PFworkjob15〜PFworkjob21	父親の職種
w1学年〜w7学年	各Waveの学年

0.2.4 欠損値コード

このデータでは，以下の特殊な値が使用されている．

コード	意味
9999	無回答・不明
8888	非該当
7777	質問なし（対象外の設問）
システム欠損	調査対象外（小学校入学前・高校卒業後）

0.2.5 多重回答の処理

多重回答の変数は以下のようにコーディングされている．

コード	意味
1	選択
0	非選択
9999	無回答・不明

0.2.6 公開データの制限

公開データでは以下の情報が含まれていない．

• 地域情報（都道府県，市区町村など）

注意

欠損値処理を行う前に，各変数で7777, 8888, 9999がどのような意味を持つかを確認すること．変数によっては実際の値として使用されている可能性がある．

0.3 本日の流れ

時間	内容
10:30-12:00	準備編：Rの操作，プロジェクト設定，データ読み込み
12:00-13:00	昼休み
13:00-15:00	基礎編：変数の変換，記述統計，可視化
15:00-15:15	休憩
15:15-17:00	実践編：クロス表分析，回帰分析，パネルデータ

0.4 表記上の注意

• パッケージ：パッケージ名パッケージ（例：tidyverseパッケージ）
• 関数：関数名()（例：mean()関数）
• 特定のパッケージの関数：パッケージ名::関数名()（例：dplyr::filter()）

本資料では，内容に応じて4種類の色分けされた枠（callout）を使用する．

ヒント

実践的なコツやテクニックを記載する．

補足

補足情報や背景説明を記載する．

注意

間違えやすい点や落とし穴を記載する．

重要

必ず理解すべき重要な概念を記載する．

1 準備編

1.1 RStudioの起動

RStudioを起動すると，4つの領域（pane）が表示される．

• 左上：Source pane：スクリプトを表示・編集する
• 左下：Console pane：結果が表示される
• 右上：Environment pane：作成されたオブジェクトを確認できる
• 右下：Output pane：ファイル，図，パッケージ，ヘルプなどが表示される

1.1.1 よく使うショートカットキー

操作	Windows	Mac
コードを実行	Ctrl + Enter	Cmd + Enter
すべてのコードを実行	Ctrl + Shift + Enter	Cmd + Shift + Enter
パイプ演算子 |>	Ctrl + Shift + M	Cmd + Shift + M
代入演算子 <-	Alt + -	Option + -
コメントアウト	Ctrl + Shift + C	Cmd + Shift + C
セクション挿入	Ctrl + Shift + R	Cmd + Shift + R
ヘルプを表示	F1（カーソル位置の関数）	F1
矩形選択	Alt + ドラッグ	Option + ドラッグ

矩形選択の活用

矩形選択を使うと，複数行の同じ位置を同時に編集できる．パネルデータで同じ変数の複数年分を作成する際に便利である．

1. まず1行だけ書く：wake_time__2015 = ds0000000080_w1c,
2. その行を7行分コピーする
3. 矩形選択で2015の部分を選択し，2016, 2017…と書き換える
4. 同様にw1の部分もw2, w3…と書き換える

1.1.2 ヘルプの使い方

関数の使い方がわからない場合は，ヘルプを参照する．

# 関数のヘルプを表示
?mean
?lm

# 演算子のヘルプ（引用符で囲む）
?"%in%"
?"["

# キーワードで検索
help.search("regression")

Output paneの「Help」タブにヘルプが表示される．また，カーソルを関数名の上に置いてF1キーを押してもヘルプが表示される．

1.2 分析の再現性

研究において，分析結果を第三者が再現できることは非常に重要である．Rを使った分析では，コンソールに直接入力するのではなく，すべての操作をスクリプトファイル（.Rや.qmd）に記録することで，再現可能な研究（Reproducible Research）を実現できる．

スクリプトに残すことの重要性

コンソールに直接入力すると記録に残らない．必ずRスクリプトに入力し，それを実行するようにする．

分析の結果自体は保存する必要はない．分析のプロセスがRスクリプトに残っていれば，結果は何度でも再現できる．

1.3 Rによる計算

まずは簡単な計算から始める．

# 足し算
1 + 1

[1] 2

# 引き算
2 - 100

[1] -98

# 掛け算
7 * 8

[1] 56

# 割り算
123456 / 3

[1] 41152

# 累乗
2^3

[1] 8

[1]は何？

結果の左に表示される[1]は，結果をベクトルで表示した時に，画面の一番左にある数字が何番目の要素であるかを示すものである．

ESCキーの使い方

計算の途中で改行してしまったり，コードに問題があって結果が表示されない場合は，コンソールに+が表示されて入力待ちの状態になる．このような場合はESCキーを押すことで>の状態に戻すことができる．何かおかしいと思ったら，まずESCで元の状態に戻す．

#（ハッシュ）でコメントを残す

#でコメントをマメに残す癖をつける．コードを書いているときは理解できていても，時間が空くと何をやっているかを忘れてしまうことがよくある．他人にコードをチェックしてもらう際にもコメントは重要である．

# --------------------------
# 2026年2月17日
# --------------------------
# 基礎集計

# 性別の度数
d |> count(sex)

RStudioの「Code」→「Insert Section…」でセクション区切りを挿入することもできる．

1.4 関数の使用

関数は関数名(引数)の形で使用する．

# 平方根
sqrt(8)

[1] 2.828427

# 自然対数
log(8)

[1] 2.079442

# 底が2の対数
log(8, base = 2)

[1] 3

# 指数関数
exp(2)

[1] 7.389056

1.5 ベクトル

複数の数値を並べるにはc()を使用する．

# ベクトルの作成
c(1, 2, 2, 3)

[1] 1 2 2 3

# 連続した値
1:10

 [1]  1  2  3  4  5  6  7  8  9 10

# 繰り返し
rep(1, times = 5)

[1] 1 1 1 1 1

1.6 数列の作成：seq()

seq()関数は等差数列を作成する．

# 0から10まで2刻み
seq(0, 10, by = 2)

[1]  0  2  4  6  8 10

# 0から10まで6個の数を等間隔で
seq(0, 10, length.out = 6)

[1]  0  2  4  6  8 10

# 1から始めて，2ずつ増やして5個
seq(from = 1, by = 2, length.out = 5)

[1] 1 3 5 7 9

seq()の主な引数

• from：開始値
• to：終了値
• by：増分
• length.out：要素数

1.7 無作為抽出：sample()

sample()関数はベクトルから無作為に要素を抽出する．

# 1から10の中から3つを無作為に抽出（非復元抽出）
sample(1:10, size = 3)

[1] 9 4 8

# 復元抽出（同じ値が複数回選ばれる可能性あり）
sample(1:10, size = 5, replace = TRUE)

[1]  6  5  1  8 10

# 文字ベクトルからの抽出
sample(c("A", "B", "C", "D"), size = 2)

[1] "C" "A"

# ベクトルの順番をシャッフル
sample(1:5)

[1] 3 4 5 2 1

再現性の確保：set.seed()

sample()は実行するたびに結果が変わる．同じ結果を再現したい場合は，直前にset.seed()で乱数の種を固定する．

set.seed(123)
sample(1:10, 3)  # 何度実行しても同じ結果

データフレームからの無作為抽出：slice_sample()

データフレームから無作為に行を抽出する場合は，dplyrのslice_sample()を使用する．

# 100行を無作為抽出
d |> slice_sample(n = 100)

# 10%を無作為抽出
d |> slice_sample(prop = 0.1)

# 復元抽出（ブートストラップなどで使用）
d |> slice_sample(n = 100, replace = TRUE)

1.8 文字列の結合：paste()とpaste0()

paste()とpaste0()は文字列を結合する．

# paste()：スペースで結合（デフォルト）
paste("Hello", "World")

[1] "Hello World"

# sep引数で区切り文字を指定
paste("2024", "01", "15", sep = "-")

[1] "2024-01-15"

# paste0()：区切りなしで結合（sep = ""と同じ）
paste0("変数", 1:3)

[1] "変数1" "変数2" "変数3"

# collapse引数でベクトルの要素を1つの文字列に
paste(c("A", "B", "C"), collapse = ", ")

[1] "A, B, C"

# sepとcollapseの組み合わせ
paste0("item_", 1:3, collapse = " + ")

[1] "item_1 + item_2 + item_3"

pasteとpaste0の使い分け

• paste()：単語をスペースで区切って結合したいとき
• paste0()：変数名やファイル名など，区切りなしで結合したいとき
• sep：複数の引数の間の区切り文字
• collapse：ベクトルの要素を1つの文字列にまとめるときの区切り文字

1.9 文字列の検索：grepl()とgrep()

文字列のパターンマッチングにはgrepl()とgrep()を使う．

# サンプルデータ
fruits <- c("apple", "banana", "cherry", "apricot", "blueberry")

# grepl()：パターンに一致するかどうか（TRUE/FALSE）
grepl("a", fruits)

[1]  TRUE  TRUE FALSE  TRUE FALSE

# grep()：パターンに一致する要素の位置（インデックス）
grep("a", fruits)

[1] 1 2 4

# パターンに一致する要素を取り出す
fruits[grepl("^a", fruits)]  # "a"で始まる要素

[1] "apple"   "apricot"

fruits[grepl("y$", fruits)]  # "y"で終わる要素

[1] "cherry"    "blueberry"

1.9.1 正規表現の基本

パターンマッチングでは正規表現（regular expression）が使える．

パターン	意味	例
^	文字列の先頭	^a：aで始まる
$	文字列の末尾	y$：yで終わる
.	任意の1文字	a.c：aとcの間に1文字
*	直前の文字が0回以上	ab*：a, ab, abb, …
+	直前の文字が1回以上	ab+：ab, abb, …
[abc]	a, b, cのいずれか	[aeiou]：母音
[^abc]	a, b, c以外	[^0-9]：数字以外
\|	または	apple\|banana

# starwarsデータで練習
library(tidyverse)

# 名前に"Sky"を含むキャラクター
starwars$name[grepl("Sky", starwars$name)]

[1] "Luke Skywalker"   "Anakin Skywalker" "Shmi Skywalker"  

# 名前が"a"で終わるキャラクター
starwars$name[grepl("a$", starwars$name)]

 [1] "Leia Organa"         "Chewbacca"           "Yoda"               
 [4] "Mon Mothma"          "Padmé Amidala"      "Sebulba"            
 [7] "Quarsh Panaka"       "Bib Fortuna"         "Ayla Secura"        
[10] "Adi Gallia"          "Mas Amedda"          "Bail Prestor Organa"
[13] "Captain Phasma"     

stringrパッケージ

tidyverseに含まれるstringrパッケージは，より直感的な文字列処理関数を提供する．

library(stringr)
str_detect(fruits, "a")      # grepl()と同じ
str_subset(fruits, "^a")     # パターンに一致する要素を返す
str_replace(fruits, "a", "X") # 最初のaをXに置換
str_replace_all(fruits, "a", "X") # すべてのaをXに置換

1.10 便利な関数（データ分析での活用例）

以下の関数は，後のデータ分析で頻繁に使用する．ここでは基本的な使い方を紹介し，実際のデータでの活用場面を示す．

1.10.1 データの大きさを確認：nrow()，ncol()，length()

データを読み込んだ後，まずサンプルサイズと変数の数を確認することが重要である．

# サンプルデータで練習
df <- data.frame(
  id = 1:5,
  score = c(80, 90, 75, 85, 95)
)

# 行数（サンプルサイズ）
nrow(df)

[1] 5

# 列数（変数の数）
ncol(df)

[1] 2

# 行数と列数を同時に取得
dim(df)

[1] 5 2

# ベクトルの長さ
length(1:10)

[1] 10

データ分析での活用場面

• nrow(d) でサンプルサイズを確認
• ncol(d) で変数の数を確認（このデータは1000以上の変数がある）
• 分析対象を絞った後に nrow() で何ケース残っているか確認

length()関数の注意点

欠損値がある場合，length()関数を有効ケース数として使用してはいけない．length()は欠損値を含む全要素数を返す．

x <- c(1, 2, NA, 4, NA)
length(x)           # 5（欠損を含む）
sum(!is.na(x))      # 3（欠損を除いた有効ケース数）

有効ケース数を求めるにはsum(!is.na(x))を使用する．頻繁に使う場合は関数を作成すると便利である．

# 欠損値を除いたサンプルサイズを求める関数
complete_obs <- function(x) sum(!is.na(x))
complete_obs(x)  # 3

nrow()も同様に欠損値を考慮しないため，特定の変数に欠損がないケース数を知りたい場合は注意が必要である．

1.10.2 重複を除く：unique()

変数がどのような値を取るか確認するのに便利である．

# 重複のあるベクトル
x <- c(1, 2, 2, 3, 3, 3, 4)

# 重複を除いた値
unique(x)

[1] 1 2 3 4

# ユニークな値の数を数える
length(unique(x))

[1] 4

データ分析での活用場面

• unique(d$w1学年) で学年変数がどのような値を取るか確認
• length(unique(d$PanelID)) でユニークな個人数を確認
• カテゴリ変数のカテゴリ数を把握するのに便利

1.10.3 並び替え：sort()

データの分布を把握するのに役立つ．

# ベクトルの並び替え
x <- c(3, 1, 4, 1, 5, 9, 2, 6)

# 昇順に並び替え
sort(x)

[1] 1 1 2 3 4 5 6 9

# 降順に並び替え
sort(x, decreasing = TRUE)

[1] 9 6 5 4 3 2 1 1

データ分析での活用場面

• sort(unique(d$birth_year)) で出生年度の範囲を確認
• 外れ値の確認に sort(d$income__2015) の最初と最後を見る
• データフレームの並び替えには dplyr::arrange() を使用

1.10.4 条件を満たす位置を取得：which()

特定の条件を満たすケースを探すのに使う．

x <- c(10, 25, 30, 15, 40, 5)

# 30以上の値がある位置（インデックス）
which(x >= 30)

[1] 3 5

# 最大値・最小値の位置
which.max(x)

[1] 5

which.min(x)

[1] 6

データ分析での活用場面

• which(d$vocab_w2_g3 > 70) で語彙力スコアが高いケースの行番号を取得
• which.max(d$income__2015) で最高収入のケースを特定
• 外れ値や異常値の位置を特定するのに便利

1.10.5 条件判定：any()とall()

データの特性を素早くチェックするのに使う．

x <- c(1, 2, 3, 4, 5)

# いずれかが条件を満たすか
any(x > 4)   # TRUE

[1] TRUE

# すべてが条件を満たすか
all(x > 0)   # TRUE

[1] TRUE

# 欠損値の確認
y <- c(1, 2, NA, 4)
any(is.na(y))   # NAがあるか？

[1] TRUE

all(!is.na(y))  # すべて非欠損か？

[1] FALSE

データ分析での活用場面

• any(is.na(d$vocab_w2_g3)) で欠損値があるか確認
• all(d$w1学年 %in% 1:12) ですべてのケースが在学中か確認
• データの前提条件を素早くチェック

1.10.6 欠損値の判定：is.na()

欠損値の処理はデータ分析で最も重要な作業の一つである．

x <- c(1, NA, 3, NA, 5)

# 欠損値かどうか
is.na(x)

[1] FALSE  TRUE FALSE  TRUE FALSE

# 欠損値の数をカウント
sum(is.na(x))

[1] 2

# 欠損値を除いたベクトル
x[!is.na(x)]

[1] 1 3 5

データ分析での活用場面

• sum(is.na(d$vocab_w2_g3)) で語彙力スコアの欠損数を確認
• mean(is.na(d$income__2015)) で欠損率を計算（0.3なら30%欠損）
• d[!is.na(d$vocab_w2_g3), ] で欠損のないケースだけ抽出

1.10.7 文字列操作：paste()での変数名作成

繰り返し処理や変数名の作成で活躍する（paste()は前述）．

# Wave1からWave7の変数名を一括作成
paste0("income__", 2015:2021)

[1] "income__2015" "income__2016" "income__2017" "income__2018" "income__2019"
[6] "income__2020" "income__2021"

# 複数の変数名パターンを作成
paste0("ds0000000121_w", 1:7, "p")

[1] "ds0000000121_w1p" "ds0000000121_w2p" "ds0000000121_w3p" "ds0000000121_w4p"
[5] "ds0000000121_w5p" "ds0000000121_w6p" "ds0000000121_w7p"

データ分析での活用場面

• 複数Waveの変数名を一括で指定
• ファイル名の自動生成
• ループ処理での変数名の動的指定

1.10.8 数値の丸め：round()

結果を見やすく表示するのに使う．

x <- 3.14159

# 四捨五入（小数点以下の桁数を指定）
round(x, digits = 2)

[1] 3.14

# 切り捨て・切り上げ
floor(3.9)   # 3

[1] 3

ceiling(3.1) # 4

[1] 4

データ分析での活用場面

• round(mean(d$vocab_w2_g3, na.rm = TRUE), 2) で平均値を小数点2桁に
• 表やグラフのラベルを見やすくする
• パーセントを整数に丸める

1.10.9 その他の計算に使う関数

# 合計
sum(1:5)

[1] 15

# 絶対値
abs(-5)

[1] 5

# 累積和
cumsum(1:5)

[1]  1  3  6 10 15

データ分析での活用場面

• sum(!is.na(d$vocab_w2_g3)) で有効ケース数を計算
• abs(d$score1 - d$score2) でスコアの差の絶対値を取得
• cumsum() で累積度数を計算

1.11 オブジェクトへの代入

Rでは、値や計算結果に名前をつけて保存できる．この「名前をつけた入れ物」をオブジェクトと呼ぶ．<-（代入演算子）を使って、右側の値を左側の名前に代入する．

# aに4を代入
a <- 4
a

[1] 4

# bにベクトルを代入
b <- c(1, 2, 3, 4, 5)
b

[1] 1 2 3 4 5

# 平均値
mean(b)

[1] 3

# 標準偏差
sd(b)

[1] 1.581139

<- と =

代入には<-を使用し，関数の引数の指定には=を使用するのが一般的である．

オブジェクトを活用するメリット

オブジェクトを活用することで，効率的な分析が可能になる：

• 再利用：一度読み込んだデータや計算結果を何度でも使える
• 可読性：mean(income)のように，何を計算しているか分かりやすくなる
• 修正が容易：データの前処理を修正しても，後続の分析コードを変更する必要がない
• 再現性：スクリプトを実行すれば，同じ結果を何度でも再現できる

1.12 パッケージのインストールと読み込み

本コースで使用するパッケージをインストールする．

# インストール（初回のみ）
install.packages("tidyverse", dependencies = TRUE)
install.packages("haven", dependencies = TRUE)
install.packages("janitor", dependencies = TRUE)
install.packages("here", dependencies = TRUE)
install.packages("broom", dependencies = TRUE)
install.packages("pacman", dependencies = TRUE)

パッケージの読み込みは毎回必要である．

# パッケージの読み込み
library(tidyverse)
library(haven)
library(janitor)
library(here)
library(broom)

# tibbleの表示オプション（出力が長くなりすぎないように）
options(
  tibble.print_max = 20,
  tibble.print_min = 10,
  tibble.width = 80
)

1.12.1 pacmanパッケージを使う方法

pacmanパッケージのp_load()関数を使うと，インストールと読み込みを同時に行える．パッケージがインストールされていなければ自動的にインストールし，読み込む．

# pacmanがインストールされていない場合は先にインストール
# install.packages("pacman")

# インストール（未インストールの場合のみ）と読み込みを同時に行う
pacman::p_load(tidyverse,
               haven,
               janitor,
               here,
               broom)

p_load()の利点

• インストール済みかどうかを気にしなくてよい
• 複数パッケージを1行で読み込める
• スクリプトの共有時に便利（相手の環境でも自動インストールされる）

::（ダブルコロン）の使い方

パッケージ名::関数名()の形式で，パッケージを読み込まずに関数を使用できる．

# library()で読み込んでから使う
library(dplyr)
filter(d, x > 0)

# ::で直接使う（library()不要）
dplyr::filter(d, x > 0)

::の利点

• パッケージを読み込まなくても関数を使える
• どのパッケージの関数かが明確になる
• 同名関数の衝突を避けられる（例：dplyr::filter() vs stats::filter()）
• スクリプトの可読性が向上する

1.13 プロジェクトの作成

RStudioでプロジェクトを作成すると，作業ディレクトリの管理が楽になる．

1. 「File」→「New Project…」を選択
2. 「New Directory」または「Existing Directory」を選択
3. フォルダを指定して「Create Project」

プロジェクトを作成すると，.Rprojファイルが作成される．このファイルをダブルクリックすれば，自動的に作業ディレクトリが設定される．

作業スペースの保存設定（推奨）

RStudioを終了する際に「Save Workspace…」と聞かれるが，保存しないことを推奨する．以下の設定を行うこと．

「Tools」→「Global Options」→「General」で：

• 「Restore .RData into workspace at startup」のチェックを外す
• 「Save workspace to .RData on exit:」を「Never」に設定

スクリプトさえ残っていれば，いつでも結果を再現できる．作業スペースを保存すると，古いオブジェクトが残って混乱の原因になることがある．

1.13.1 hereパッケージによるパス管理

hereパッケージを使うと，プロジェクトのルートディレクトリを基準にファイルパスを指定できる．

library(here)

# プロジェクトのルートディレクトリを確認
here()

[1] "/Users/sf/statsemi_R"

# ファイルパスの指定例
here("data", "raw", "1571", "1571.dta")

[1] "/Users/sf/statsemi_R/data/raw/1571/1571.dta"

here()の利点

• OSによるパス区切り文字の違い（/ vs \）を気にしなくてよい
• setwd()を使う必要がない
• 他の人とコードを共有しても，各自の環境で動作する

1.14 フォルダ構成

以下のようなフォルダ構成を推奨する．

benesse_seminar/
├── benesse_seminar.Rproj
├── data/
│   ├── raw/          # 元データ
│   │   ├── 1571/     # メインのパネルデータ
│   │   ├── 1577/     # 2016年度語彙力・読解力調査
│   │   └── 1578/     # 2019年度語彙力・読解力調査
│   └── processed/    # 加工済みデータ
├── scripts/          # Rスクリプト
├── figures/          # 図
└── tables/           # 表

# フォルダの作成
library(fs)
dir_create(here(c("data/raw",
                  "data/processed",
                  "scripts",
                  "figures",
                  "tables")))

1.15 データの読み込み

1.15.1 メインデータ（1571）の読み込み

# Stataファイル（.dta）で読み込む場合
d_1571 <- read_dta(here("data", "raw", "1571", "1571.dta"))

# CSVファイルで読み込む場合
# d <- read_csv(here("data", "raw", "1571", "1571.csv"))

# SPSSファイル（.sav）で読み込む場合
# d <- read_sav(here("data", "raw", "1571", "1571.sav"))

# データの確認
dim(d_1571)  # 行数と列数

[1] 29846  4028

Stataファイルのラベル属性について

haven::read_dta()で読み込んだStataファイルの変数は，値ラベル（value labels）の情報をhaven_labelledクラスとして保持している．この属性があるとdplyr::case_match()などで予期しない動作をすることがあるため，値をリコードする際はhaven::zap_labels()でラベル属性を除去する．

すべての変数から一括でラベルを除去したい場合は，読み込み直後に以下を実行する：

d_1571 <- d_1571 |> mutate(across(everything(), zap_labels))

一方，値ラベルをそのままfactor型のカテゴリ名として使いたい場合はhaven::as_factor()が便利である：

# 値ラベルをfactorのlevelsに変換
d <- d |> mutate(var_f = as_factor(var))

as_factor()を使うと値ラベルがそのままカテゴリ名になるため，case_match()で手動マッピングする手間が省ける．ただし「非該当」「無回答」などもカテゴリに含まれるため，na_if()やfct_drop()で適宜処理する必要がある．

1.15.2 語彙力・読解力調査データの読み込み

# 2016年度（Wave2時点）語彙力・読解力調査
d_1577 <- read_dta(here("data", "raw", "1577", "1577.dta"))
# 変数名を確認
# names(d_1577)
# 語彙力・読解力の偏差値スコアを抽出してリネーム
d_1577 <- d_1577 |>
  select(PanelID,
         vocab_w2_g3 = contains("語彙w2") & contains("小３偏差値"),
         vocab_w2_g9 = contains("読解w2") & contains("中３偏差値"))

# 2019年度（Wave5時点）語彙力・読解力調査
d_1578 <- read_dta(here("data", "raw", "1578", "1578.dta"))
d_1578 <- d_1578 |>
  select(PanelID,
         vocab_w5_g3 = contains("語彙w5") & contains("小３偏差値"),
         vocab_w5_g9 = contains("読解w5") & contains("中３偏差値"))

1.15.2.1 語彙力・読解力データの変数

変数名	内容
PanelID	パネルID（メインデータとの結合キー）
語彙wX_学年	語彙力調査に回答した時点の学年
語彙wX_回答フラグ	語彙力調査に協力した者を1
語彙wX_IRTスコア	語彙力調査のIRTスコア
語彙wX_小３偏差値スコア	IRTスコアの小3基準偏差値（小3の平均が50）
読解wX_学年	読解力調査に回答した時点の学年
読解wX_回答フラグ	読解力調査に協力した者を1
読解wX_IRTスコア	読解力調査のIRTスコア
読解wX_中３偏差値スコア	IRTスコアの中3基準偏差値（中3の平均が50）

1.15.3 データの結合

メインデータと語彙力・読解力データをPanelIDで結合する．

# メインデータに語彙力・読解力データを結合
d <- d_1571 |> 
  left_join(d_1577, by = "PanelID") |> 
  left_join(d_1578, by = "PanelID")

left_join()について

left_join()は左側のデータ（ここではメインデータd）のすべての行を保持し，右側のデータ（語彙力データ）から一致する行の情報を追加する．一致しない場合はNAが入る．

1.15.4 欠損値の処理

# 特殊コードをNAに変換する関数
recode_missing <- function(x) {
  if_else(x %in% c(7777, 8888, 9999), NA, x)
}

# 特定の変数に適用
d <- d |> 
  mutate(
    across(
      c(contains("子ども性別"), starts_with("PFED"), starts_with("PMED")),
      recode_missing
    )
  )

1.15.5 変数の確認

# 変数名の一覧（最初の50個）
names(d) |> head(50)

 [1] "PanelID"      "w1回答フラグ" "w2回答フラグ" "w3回答フラグ" "w4回答フラグ"
 [6] "w5回答フラグ" "w6回答フラグ" "w7回答フラグ" "w1学年"       "w2学年"      
[11] "w3学年"       "w4学年"       "w5学年"       "w6学年"       "w7学年"      
[16] "w1子ども性別" "w2子ども性別" "w3子ども性別" "w4子ども性別" "w5子ども性別"
[21] "w6子ども性別" "w7子ども性別" "PFED15"       "PFED16"       "PFED17"      
[26] "PFED18"       "PFED19"       "PFED20"       "PFED21"       "PMED15"      
[31] "PMED16"       "PMED17"       "PMED18"       "PMED19"       "PMED20"      
[36] "PMED21"       "PFworkhome16" "PFworkhome17" "PFworkhome18" "PFworkhome19"
[41] "PFworkhome20" "PFworkhome21" "PFworkjob15"  "PFworkjob16"  "PFworkjob17" 
[46] "PFworkjob18"  "PFworkjob19"  "PFworkjob20"  "PFworkjob21"  "PFworktime15"

# 特定の変数の分布を確認
d |> count(w1学年)

# A tibble: 14 × 2
   w1学年                          n
   <dbl+lbl>                   <int>
 1    1 [小１]                  1837
 2    2 [小２]                  1577
 3    3 [小３]                  1627
 4    4 [小４]                  1465
 5    5 [小５]                  1415
 6    6 [小６]                  1420
 7    7 [中１]                  1475
 8    8 [中２]                  1539
 9    9 [中３]                  1506
10   10 [高１]                  1421
11   11 [高２]                  1420
12   12 [高３]                  1382
13 8888 [調査対象外（就学前）]  5541
14   NA                         6221

d |> count(w1子ども性別)

# A tibble: 3 × 2
  w1子ども性別     n
  <dbl+lbl>    <int>
1  1 [男子]     8300
2  2 [女子]     8307
3 NA           13239

2 基礎編

以降の分析では，準備編で読み込んだデータdを使用する．dは「子どもの生活と学びに関する親子調査」のWave 1からWave 5のデータを結合したものである．

2.1 パイプ演算子

パイプ演算子|>を使うと，処理の流れが分かりやすくなる．

# 通常の書き方
round(mean(d$vocab_w2_g3, 
           na.rm = TRUE), digits = 1)

[1] 57.4

# パイプを使った書き方
d$vocab_w2_g3 |> 
  mean(na.rm = TRUE) |> 
  round(digits = 1)

[1] 57.4

パイプの左側の結果が，右側の関数の第1引数に渡される．

2.2 データの確認

データを読み込んだら、まず全体の概要を確認する。

# glimpse()でデータの構造を確認（変数が多いため出力は省略）
glimpse(d)

# 先頭6行
head(d)

# A tibble: 6 × 4,032
  PanelID w1回答フラグ   w2回答フラグ     w3回答フラグ w4回答フラグ w5回答フラグ
    <dbl> <dbl+lbl>      <dbl+lbl>        <dbl+lbl>    <dbl+lbl>    <dbl+lbl>   
1 1600001 6 [調査対象外] 1 [親子とも回答（小１～３含… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小…
2 1600002 6 [調査対象外] 1 [親子とも回答（小１～３含… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小…
3 1600003 6 [調査対象外] 1 [親子とも回答（小１～３含… 1 [親子とも回答（小… 1 [親子とも回答（小… 4 [回答なし]   …
4 1600004 6 [調査対象外] 1 [親子とも回答（小１～３含… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小…
5 1600005 6 [調査対象外] 1 [親子とも回答（小１～３含… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小…
6 1600006 6 [調査対象外] 1 [親子とも回答（小１～３含… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小…
# ℹ 4,026 more variables: w6回答フラグ <dbl+lbl>, w7回答フラグ <dbl+lbl>,
#   w1学年 <dbl+lbl>, w2学年 <dbl+lbl>, w3学年 <dbl+lbl>, w4学年 <dbl+lbl>,
#   w5学年 <dbl+lbl>, w6学年 <dbl+lbl>, w7学年 <dbl+lbl>,
#   w1子ども性別 <dbl+lbl>, w2子ども性別 <dbl+lbl>, w3子ども性別 <dbl+lbl>,
#   w4子ども性別 <dbl+lbl>, w5子ども性別 <dbl+lbl>, w6子ども性別 <dbl+lbl>,
#   w7子ども性別 <dbl+lbl>, PFED15 <dbl+lbl>, PFED16 <dbl+lbl>,
#   PFED17 <dbl+lbl>, PFED18 <dbl+lbl>, PFED19 <dbl+lbl>, PFED20 <dbl+lbl>, …

# 要約統計量（一部の変数のみ）
d |> select(PanelID, w1子ども性別, vocab_w2_g3, vocab_w2_g9) |> summary()

    PanelID         w1子ども性別    vocab_w2_g3      vocab_w2_g9   
 Min.   :1400001   Min.   :1.0     Min.   : 27.37   Min.   :16.47  
 1st Qu.:1507454   1st Qu.:1.0     1st Qu.: 49.20   1st Qu.:44.39  
 Median :1516138   Median :2.0     Median : 56.51   Median :50.15  
 Mean   :1640509   Mean   :1.5     Mean   : 57.41   Mean   :51.73  
 3rd Qu.:1800397   3rd Qu.:2.0     3rd Qu.: 64.12   3rd Qu.:58.00  
 Max.   :9999994   Max.   :2.0     Max.   :105.23   Max.   :91.84  
                   NA's   :13239   NA's   :26213    NA's   :28730  

2.2.1 識別変数の重複チェック

パネルデータでは識別変数（ID）が一意であることを確認する。重複があると、データの結合や集計で問題が生じる。

# 識別変数の重複チェック
sum(duplicated(d$PanelID))  # 0なら重複なし

[1] 0

# 重複があった場合、該当レコードを確認
# d |> filter(duplicated(PanelID) | duplicated(PanelID, fromLast = TRUE))

データチェックの考え方

すべての変数を事前にチェックするのは現実的ではない。実務では：

1. 読み込み直後：行数・列数、識別変数の重複を確認
2. 分析に使う変数を選んだ後：その変数の分布・欠損・値ラベルを確認（count()、summary()）
3. 変数作成後：変換前後の対応を確認（クロス表）

「分析の大半はデータの前処理と変数作成である」という意識を持ち、各段階で確認を怠らないことが重要。

2.3 変数へのアクセス

データフレーム内の変数には$でアクセスする．

# 語彙力スコアという変数
d$vocab_w2_g3 |> head(50)

 [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
[26] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA

# 語彙力スコアの平均値
mean(d$vocab_w2_g3, na.rm = TRUE)

[1] 57.41331

2.4 dplyrによるデータ操作

dplyrパッケージの主要な関数を紹介する．

2.4.1 count()：度数を数える

# 性別の度数
d |> count(w1子ども性別)

# A tibble: 3 × 2
  w1子ども性別     n
  <dbl+lbl>    <int>
1  1 [男子]     8300
2  2 [女子]     8307
3 NA           13239

# 学年の度数
d |> count(w1学年)

# A tibble: 14 × 2
   w1学年                          n
   <dbl+lbl>                   <int>
 1    1 [小１]                  1837
 2    2 [小２]                  1577
 3    3 [小３]                  1627
 4    4 [小４]                  1465
 5    5 [小５]                  1415
 6    6 [小６]                  1420
 7    7 [中１]                  1475
 8    8 [中２]                  1539
 9    9 [中３]                  1506
10   10 [高１]                  1421
11   11 [高２]                  1420
12   12 [高３]                  1382
13 8888 [調査対象外（就学前）]  5541
14   NA                         6221

# クロス集計
d |> count(w1子ども性別, w1学年)

# A tibble: 38 × 3
   w1子ども性別 w1学年        n
   <dbl+lbl>    <dbl+lbl> <int>
 1 1 [男子]      1 [小１]   924
 2 1 [男子]      2 [小２]   716
 3 1 [男子]      3 [小３]   725
 4 1 [男子]      4 [小４]   647
 5 1 [男子]      5 [小５]   635
 6 1 [男子]      6 [小６]   655
 7 1 [男子]      7 [中１]   653
 8 1 [男子]      8 [中２]   673
 9 1 [男子]      9 [中３]   694
10 1 [男子]     10 [高１]   637
# ℹ 28 more rows

2.4.2 filter()：行の抽出

# 男子（性別=1）のみ
d |> filter(w1子ども性別 == 1)

# A tibble: 8,300 × 4,032
   PanelID w1回答フラグ      w2回答フラグ w3回答フラグ w4回答フラグ w5回答フラグ
     <dbl> <dbl+lbl>         <dbl+lbl>    <dbl+lbl>    <dbl+lbl>    <dbl+lbl>   
 1 9999991 1 [親子とも回答（小１～３含）… 4 [回答なし]   … 6 [調査対象外]… … 6 [調査対象外]… … 6 [調査対象外]… …
 2 9999992 1 [親子とも回答（小１～３含）… 4 [回答なし]   … 6 [調査対象外]… … 6 [調査対象外]… … 6 [調査対象外]… …
 3 9999994 1 [親子とも回答（小１～３含）… 4 [回答なし]   … 6 [調査対象外]… … 6 [調査対象外]… … 6 [調査対象外]… …
 4 1517967 1 [親子とも回答（小１～３含）… 1 [親子とも回答（小… 1 [親子とも回答（小… 4 [回答なし]   … 4 [回答なし]   …
 5 1515102 1 [親子とも回答（小１～３含）… 1 [親子とも回答（小… 4 [回答なし]   … 4 [回答なし]   … 4 [回答なし]   …
 6 1514527 1 [親子とも回答（小１～３含）… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小…
 7 1516893 1 [親子とも回答（小１～３含）… 4 [回答なし]   … 4 [回答なし]   … 4 [回答なし]   … 4 [回答なし]   …
 8 1514307 1 [親子とも回答（小１～３含）… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小…
 9 1517985 1 [親子とも回答（小１～３含）… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小…
10 1517952 1 [親子とも回答（小１～３含）… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小…
# ℹ 8,290 more rows
# ℹ 4,026 more variables: w6回答フラグ <dbl+lbl>, w7回答フラグ <dbl+lbl>,
#   w1学年 <dbl+lbl>, w2学年 <dbl+lbl>, w3学年 <dbl+lbl>, w4学年 <dbl+lbl>,
#   w5学年 <dbl+lbl>, w6学年 <dbl+lbl>, w7学年 <dbl+lbl>,
#   w1子ども性別 <dbl+lbl>, w2子ども性別 <dbl+lbl>, w3子ども性別 <dbl+lbl>,
#   w4子ども性別 <dbl+lbl>, w5子ども性別 <dbl+lbl>, w6子ども性別 <dbl+lbl>,
#   w7子ども性別 <dbl+lbl>, PFED15 <dbl+lbl>, PFED16 <dbl+lbl>, …

# 中学生（学年7-9）のみ
d |> filter(w1学年 %in% 7:9)

# A tibble: 4,520 × 4,032
   PanelID w1回答フラグ      w2回答フラグ w3回答フラグ w4回答フラグ w5回答フラグ
     <dbl> <dbl+lbl>         <dbl+lbl>    <dbl+lbl>    <dbl+lbl>    <dbl+lbl>   
 1 9999994 1 [親子とも回答（小１～３含）… 4 [回答なし]   … 6 [調査対象外]… … 6 [調査対象外]… … 6 [調査対象外]… …
 2 1505411 1 [親子とも回答（小１～３含）… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小…
 3 1505503 1 [親子とも回答（小１～３含）… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小…
 4 1511683 4 [回答なし]        … 4 [回答なし]   … 1 [親子とも回答（小… 4 [回答なし]   … 4 [回答なし]   …
 5 1509672 5 [分析対象外]       … 2 [親のみ回答]… … 1 [親子とも回答（小… 1 [親子とも回答（小… 5 [分析対象外]… …
 6 1509905 4 [回答なし]        … 1 [親子とも回答（小… 1 [親子とも回答（小… 4 [回答なし]   … 5 [分析対象外]… …
 7 1507527 1 [親子とも回答（小１～３含）… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小… 4 [回答なし]   …
 8 1509597 3 [子どものみ回答]…    … 1 [親子とも回答（小… 4 [回答なし]   … 1 [親子とも回答（小… 4 [回答なし]   …
 9 1510531 4 [回答なし]        … 1 [親子とも回答（小… 1 [親子とも回答（小… 6 [調査対象外]… … 4 [回答なし]   …
10 1505614 1 [親子とも回答（小１～３含）… 1 [親子とも回答（小… 1 [親子とも回答（小… 6 [調査対象外]… … 1 [親子とも回答（小…
# ℹ 4,510 more rows
# ℹ 4,026 more variables: w6回答フラグ <dbl+lbl>, w7回答フラグ <dbl+lbl>,
#   w1学年 <dbl+lbl>, w2学年 <dbl+lbl>, w3学年 <dbl+lbl>, w4学年 <dbl+lbl>,
#   w5学年 <dbl+lbl>, w6学年 <dbl+lbl>, w7学年 <dbl+lbl>,
#   w1子ども性別 <dbl+lbl>, w2子ども性別 <dbl+lbl>, w3子ども性別 <dbl+lbl>,
#   w4子ども性別 <dbl+lbl>, w5子ども性別 <dbl+lbl>, w6子ども性別 <dbl+lbl>,
#   w7子ども性別 <dbl+lbl>, PFED15 <dbl+lbl>, PFED16 <dbl+lbl>, …

# 語彙力スコアが60以上
d |> filter(vocab_w2_g3 >= 60)

# A tibble: 1,357 × 4,032
   PanelID w1回答フラグ      w2回答フラグ w3回答フラグ w4回答フラグ w5回答フラグ
     <dbl> <dbl+lbl>         <dbl+lbl>    <dbl+lbl>    <dbl+lbl>    <dbl+lbl>   
 1 1516906 1 [親子とも回答（小１～３含）… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小…
 2 1516903 1 [親子とも回答（小１～３含）… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小…
 3 1516909 1 [親子とも回答（小１～３含）… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小…
 4 1516951 4 [回答なし]        … 1 [親子とも回答（小… 4 [回答なし]   … 1 [親子とも回答（小… 1 [親子とも回答（小…
 5 1515950 1 [親子とも回答（小１～３含）… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小…
 6 1516669 4 [回答なし]        … 4 [回答なし]   … 4 [回答なし]   … 4 [回答なし]   … 4 [回答なし]   …
 7 1516932 1 [親子とも回答（小１～３含）… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小…
 8 1516927 1 [親子とも回答（小１～３含）… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小…
 9 1516923 1 [親子とも回答（小１～３含）… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小…
10 1516594 1 [親子とも回答（小１～３含）… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小… 1 [親子とも回答（小…
# ℹ 1,347 more rows
# ℹ 4,026 more variables: w6回答フラグ <dbl+lbl>, w7回答フラグ <dbl+lbl>,
#   w1学年 <dbl+lbl>, w2学年 <dbl+lbl>, w3学年 <dbl+lbl>, w4学年 <dbl+lbl>,
#   w5学年 <dbl+lbl>, w6学年 <dbl+lbl>, w7学年 <dbl+lbl>,
#   w1子ども性別 <dbl+lbl>, w2子ども性別 <dbl+lbl>, w3子ども性別 <dbl+lbl>,
#   w4子ども性別 <dbl+lbl>, w5子ども性別 <dbl+lbl>, w6子ども性別 <dbl+lbl>,
#   w7子ども性別 <dbl+lbl>, PFED15 <dbl+lbl>, PFED16 <dbl+lbl>, …

2.4.3 select()：列の選択

# 特定の列を選択
d |> select(PanelID, w1子ども性別, vocab_w2_g3)

# A tibble: 29,846 × 3
   PanelID w1子ども性別 vocab_w2_g3
     <dbl> <dbl+lbl>          <dbl>
 1 1600001 NA                    NA
 2 1600002 NA                    NA
 3 1600003 NA                    NA
 4 1600004 NA                    NA
 5 1600005 NA                    NA
 6 1600006 NA                    NA
 7 1600007 NA                    NA
 8 1600008 NA                    NA
 9 1600009 NA                    NA
10 1600010 NA                    NA
# ℹ 29,836 more rows

# 特定のパターンで始まる列を選択
d |> select(PanelID, starts_with("w1"))

# A tibble: 29,846 × 4
   PanelID w1回答フラグ   w1学年                      w1子ども性別
     <dbl> <dbl+lbl>      <dbl+lbl>                   <dbl+lbl>   
 1 1600001 6 [調査対象外] 8888 [調査対象外（就学前）] NA          
 2 1600002 6 [調査対象外] 8888 [調査対象外（就学前）] NA          
 3 1600003 6 [調査対象外] 8888 [調査対象外（就学前）] NA          
 4 1600004 6 [調査対象外] 8888 [調査対象外（就学前）] NA          
 5 1600005 6 [調査対象外] 8888 [調査対象外（就学前）] NA          
 6 1600006 6 [調査対象外] 8888 [調査対象外（就学前）] NA          
 7 1600007 6 [調査対象外] 8888 [調査対象外（就学前）] NA          
 8 1600008 6 [調査対象外] 8888 [調査対象外（就学前）] NA          
 9 1600009 6 [調査対象外] 8888 [調査対象外（就学前）] NA          
10 1600010 6 [調査対象外] 8888 [調査対象外（就学前）] NA          
# ℹ 29,836 more rows

2.4.4 mutate()：新しい変数の作成

# 新しい変数を作成
d <- d |> 
  mutate(
    vocab_z = (vocab_w2_g3 - 
                 mean(vocab_w2_g3, na.rm = TRUE)) / 
      sd(vocab_w2_g3, na.rm = TRUE)
  )

d |> select(PanelID, vocab_w2_g3, vocab_z) |> 
  drop_na(vocab_z)  # vocab_zに欠損がないケースのみを表示

# A tibble: 3,633 × 3
   PanelID vocab_w2_g3 vocab_z
     <dbl>       <dbl>   <dbl>
 1 1516123        44.7  -1.11 
 2 1516904        55.9  -0.129
 3 1516906        60.6   0.279
 4 1516900        41.6  -1.38 
 5 1516238        40.2  -1.50 
 6 1516498        42.3  -1.32 
 7 1516903        83.1   2.25 
 8 1516758        42.1  -1.34 
 9 1516909        64.3   0.600
10 1516136        51.0  -0.564
# ℹ 3,623 more rows

2.4.5 変数のリコード：case_match()とcase_when()

case_match vs case_when

• case_match()：特定の値を別の値に置き換える（シンプル）
• case_when()：条件式を使った柔軟な分岐が可能

# case_match：値の直接的な置換（性別）
d <- d |> 
  mutate(
    sex = coalesce(w1子ども性別,
                   w2子ども性別,
                   w3子ども性別,
                   w4子ども性別,
                   w5子ども性別,
                   w6子ども性別,
                   w7子ども性別))
d <- d |> 
  mutate(
    sex = case_match(
      sex,
      1 ~ "男子",
      2 ~ "女子",
      .default = NA
    )
  )
d |> count(sex)

# A tibble: 3 × 2
  sex       n
  <chr> <int>
1 女子  14607
2 男子  14684
3 <NA>    555

# case_match：父学歴の再カテゴリ化
d <- d |> 
  mutate(
    PFED = coalesce(zap_labels(PFED15), 
                    zap_labels(PFED16), 
                    zap_labels(PFED17), 
                    zap_labels(PFED18), 
                    zap_labels(PFED19), 
                    zap_labels(PFED20),
                    zap_labels(PFED21))
  )

d <- d |> 
  mutate(f_edu = case_match(
    PFED,
      c(1, 2) ~ "中学・高校",
      c(3, 4) ~ "専門・短大",
      c(5, 6) ~ "大学・院",
      .default = NA
    )
  )

# factorあるいはfct_relevelで水準に順位をつけたほうがよい
d <- d |> 
  mutate(f_edu = fct_relevel(f_edu, "中学・高校", "専門・短大"))


# 母親の学歴
d <- d |> 
  mutate(
    PMED = coalesce(zap_labels(PMED15),
                    zap_labels(PMED16), 
                    zap_labels(PMED17), 
                    zap_labels(PMED18), 
                    zap_labels(PMED19), 
                    zap_labels(PMED20), 
                    zap_labels(PMED21)))
d <- d |> 
  mutate(m_edu = case_match(
      PMED,
      c(1, 2) ~ "中学・高校",
      c(3, 4) ~ "専門・短大",
      c(5, 6) ~ "大学・院",
      .default = NA
    )
  )

# factorあるいはfct_relevelで水準に順位をつけたほうがよい
d <- d |>
  mutate(m_edu = fct_relevel(m_edu, "中学・高校", "専門・短大"))


# 父親の職種（W1の値を使用）
d <- d |>
  mutate(
    f_occ = case_match(
      zap_labels(PFworkjob15),
      1 ~ "事務・営業",
      2 ~ "販売・サービス",
      3 ~ "技能・労務",
      4 ~ "保安",
      5 ~ "運輸",
      6 ~ "農林漁業",
      7 ~ "専門・技術",
      8 ~ "管理職",
      c(9, 10) ~ NA,
      .default = NA
    ),
    f_occ = factor(f_occ)
  )

# 変換前後の対応を確認
d |> count(PFworkjob15, f_occ)

# A tibble: 12 × 3
   PFworkjob15                                                       f_occ     n
   <dbl+lbl>                                                         <fct> <int>
 1    1 [事務・営業系の職業（事務員、営業社員、銀行員、集金人など）]               … 事務・営…  3724
 2    2 [販売・サービス系の職業（店主、店員、外交員、美容師、清掃、ヘルパーなど）]…… 販売・サ…  1577
 3    3 [技能・労務・作業系の職業（工場労働者、建設作業者など）]                   … 技能・労…  2663
 4    4 [保安系の職業（自衛官、警察官、ガードマンなど）]                           … 保安    274
 5    5 [運輸系の職業（トラック・タクシー運転手、郵便配達、宅配など）]…            … 運輸    697
 6    6 [農林漁業職（植木職、造園業を含む）]                                       … 農林漁業…   147
 7    7 [専門・技術系の職業（医師、弁護士、教員、エンジニアなど）]                 … 専門・技…  3324
 8    8 [管理的職業（課長相当以上の管理職、議員など）]                             … 管理職……  1501
 9    9 [その他]                                                      … <NA>    832
10   10 [わからない]                                                    … <NA>    114
11 9999 [無回答・不明]                                                   … <NA>    217
12   NA                                                            … <NA>  14776

# 母親の職種（W1の値を使用）
d <- d |>
  mutate(
    m_occ = case_match(
      zap_labels(PMworkjob15),
      1 ~ "事務・営業",
      2 ~ "販売・サービス",
      3 ~ "技能・労務",
      4 ~ "保安",
      5 ~ "運輸",
      6 ~ "農林漁業",
      7 ~ "専門・技術",
      8 ~ "管理職",
      c(9, 10) ~ NA,
      .default = NA
    ),
    m_occ = factor(m_occ)
  )

# 変換前後の対応を確認
d |> count(PMworkjob15, m_occ)

# A tibble: 12 × 3
   PMworkjob15                                                       m_occ     n
   <dbl+lbl>                                                         <fct> <int>
 1    1 [事務・営業系の職業（事務員、営業社員、銀行員、集金人など）]               … 事務・営…  3416
 2    2 [販売・サービス系の職業（店主、店員、外交員、美容師、清掃、ヘルパーなど）]…… 販売・サ…  3195
 3    3 [技能・労務・作業系の職業（工場労働者、建設作業者など）]                   … 技能・労…   702
 4    4 [保安系の職業（自衛官、警察官、ガードマンなど）]                           … 保安      7
 5    5 [運輸系の職業（トラック・タクシー運転手、郵便配達、宅配など）]…            … 運輸    119
 6    6 [農林漁業職（植木職、造園業を含む）]                                       … 農林漁業…   101
 7    7 [専門・技術系の職業（医師、弁護士、教員、エンジニアなど）]                 … 専門・技…  2596
 8    8 [管理的職業（課長相当以上の管理職、議員など）]                             … 管理職……    77
 9    9 [その他]                                                      … <NA>   1299
10   10 [わからない]                                                    … <NA>      5
11 9999 [無回答・不明]                                                   … <NA>    125
12   NA                                                            … <NA>  18204

# 父親の就業形態（W1の値を使用）
d <- d |>
  mutate(
    f_emp = case_match(
      zap_labels(PFworktype15),
      1 ~ "正規雇用",
      c(2, 3, 4) ~ "非正規雇用",
      5 ~ "自営業",
      c(6, 8) ~ NA,
      7 ~ "無職",
      .default = NA
    ),
    # 非該当（8888）は父親なしとして別カテゴリに
    f_emp = if_else(zap_labels(PFworktype15) == 8888, "父親なし", f_emp),
    f_emp = factor(f_emp)
  )

# 変換前後の対応を確認
d |> count(PFworktype15, f_emp)

# A tibble: 10 × 3
   PFworktype15                      f_emp          n
   <dbl+lbl>                         <fct>      <int>
 1    1 [正社員・正職員]             正規雇用   12827
 2    2 [パート・アルバイト]         非正規雇用   107
 3    3 [契約社員・嘱託]             非正規雇用   224
 4    4 [派遣社員]                   非正規雇用    46
 5    5 [自営業（家族従業者を含む）] 自営業      1721
 6    6 [その他]                     <NA>         130
 7    7 [無職（専業主婦など）]       無職          98
 8    8 [わからない]                 <NA>          15
 9 9999 [無回答・不明]               <NA>         682
10   NA                              <NA>       13996

# 母親の就業形態（W1の値を使用）
d <- d |>
  mutate(
    m_emp = case_match(
      zap_labels(PMworktype15),
      1 ~ "正規雇用",
      c(2, 3, 4) ~ "非正規雇用",
      5 ~ "自営業",
      c(6, 8) ~ NA,
      7 ~ "無職・専業主婦",
      .default = NA
    ),
    # 非該当（8888）は母親なしとして別カテゴリに
    m_emp = if_else(zap_labels(PMworktype15) == 8888, "母親なし", m_emp),
    m_emp = factor(m_emp)
  )

# 変換前後の対応を確認
d |> count(PMworktype15, m_emp)

# A tibble: 10 × 3
   PMworktype15                      m_emp              n
   <dbl+lbl>                         <fct>          <int>
 1    1 [正社員・正職員]             正規雇用        2859
 2    2 [パート・アルバイト]         非正規雇用      6775
 3    3 [契約社員・嘱託]             非正規雇用       722
 4    4 [派遣社員]                   非正規雇用       218
 5    5 [自営業（家族従業者を含む）] 自営業           859
 6    6 [その他]                     <NA>             208
 7    7 [無職（専業主婦など）]       無職・専業主婦  4474
 8    8 [わからない]                 <NA>               1
 9 9999 [無回答・不明]               <NA>             576
10   NA                              <NA>           13154

# 自宅の本の冊数（W2で測定）
d <- d |>
  mutate(
    # カテゴリ変数
    n_books = case_match(
      zap_labels(ds0000000609_w2p),
      1 ~ "9冊以下",
      2 ~ "10〜29冊",
      3 ~ "30〜49冊",
      4 ~ "50〜99冊",
      5 ~ "100〜199冊",
      6 ~ "200〜499冊",
      7 ~ "500冊以上",
      .default = NA
    ),
    n_books = fct_relevel(n_books,
                          "9冊以下", "10〜29冊", "30〜49冊",
                          "50〜99冊", "100〜199冊", "200〜499冊"),
    # 連続変数（各カテゴリの中間値を代入）
    n_books_num = case_match(
      zap_labels(ds0000000609_w2p),
      1 ~ 5,
      2 ~ 20,
      3 ~ 40,
      4 ~ 75,
      5 ~ 150,
      6 ~ 350,
      7 ~ 500,
      .default = NA
    )
  )

# 変換前後の対応を確認
d |> count(ds0000000609_w2p, n_books, n_books_num)

# A tibble: 9 × 4
  ds0000000609_w2p    n_books    n_books_num     n
  <dbl+lbl>           <fct>            <dbl> <int>
1    1 [9冊以下]      9冊以下              5  1206
2    2 [10～29冊]     10〜29冊            20  3748
3    3 [30～49冊]     30〜49冊            40  2971
4    4 [50～99冊]     50〜99冊            75  3315
5    5 [100～199冊]   100〜199冊         150  2443
6    6 [200～499冊]   200〜499冊         350  1503
7    7 [500冊以上]    500冊以上          500   568
8 9999 [無回答・不明] <NA>                NA   259
9   NA                <NA>                NA 13833

カテゴリ変数の連続変数化について

本の冊数のようにカテゴリで測定された変数を連続変数として扱う場合，以下の点に注意が必要である：

• 各カテゴリの中間値（代表値）を代入しているが，実際の分布は不明
• 「500冊以上」のような開放カテゴリは特に不確実性が高い
• 線形の関係を仮定することの妥当性は検討が必要

カテゴリ変数（n_books）をダミー変数として使えば各カテゴリの効果を正確に推定できる．連続変数（n_books_num）はsplines::ns()で非線形性をモデル化する際に使用できる．

社会経済的地位（SES）の指標

社会経済的地位を測定する代表的な変数として以下がある：

• 収入：世帯年収（income__2015〜income__2021）
• 学歴：父親・母親の最終学歴（f_edu, m_edu）
• 職業：父親・母親の職種（f_occ, m_occ）
• 就業形態：父親・母親の就業形態（f_emp, m_emp）
• 文化資本：自宅の本の冊数（n_books, n_books_num）

これらは相互に関連しており，どれを用いるかは分析の目的による．

2.4.6 職業・就業形態の組み合わせカテゴリ（EGP分類を参考に）

職業と就業形態を組み合わせた社会階層カテゴリを作成する。EGP（Erikson-Goldthorpe-Portocarero）分類を参考に、利用可能な変数から簡易的な階層分類を試みる。

# 父親の職業階層（EGP分類を参考にした簡易版）
d <- d |>
  mutate(
    f_class = case_when(
      # 専門・管理職（正規雇用または自営）
      f_occ == "管理職" & f_emp %in% c("正規雇用", "自営業") ~ "専門・管理職",
      f_occ == "専門・技術" & f_emp %in% c("正規雇用", "自営業") ~ "専門・管理職",
      # 事務・販売職（正規雇用）
      f_occ == "事務・営業" & f_emp == "正規雇用" ~ "事務・販売職",
      f_occ == "販売・サービス" & f_emp == "正規雇用" ~ "事務・販売職",
      # 自営業（専門・管理職以外の自営、農林漁業の自営を含む）
      f_emp == "自営業" ~ "自営業",
      f_occ == "農林漁業" & f_emp == "自営業" ~ "自営業",
      # マニュアル（技能・運輸・保安で正規雇用、農林漁業の正規雇用を含む）
      f_occ %in% c("技能・労務", "運輸", "保安") & f_emp == "正規雇用" ~ "マニュアル",
      f_occ == "農林漁業" & f_emp == "正規雇用" ~ "マニュアル",
      # 非正規雇用（農林漁業の非正規を含む）
      f_emp == "非正規雇用" ~ "非正規雇用",
      # 無職
      f_emp == "無職" ~ "無職",
      # 父親なし
      f_emp == "父親なし" ~ "父親なし",
      .default = NA
    ),
    f_class = factor(f_class, levels = c(
      "専門・管理職", "事務・販売職", "自営業",
      "マニュアル", "非正規雇用", "無職", "父親なし"
    ))
  )

# 分布を確認
d |> count(f_class)

# A tibble: 7 × 2
  f_class          n
  <fct>        <int>
1 専門・管理職  4734
2 事務・販売職  4731
3 自営業        1218
4 マニュアル    3060
5 非正規雇用     377
6 無職            98
7 <NA>         15628

# 職業×就業形態の組み合わせを確認
d |> count(f_occ, f_emp, f_class) |> print(n = 50)

# A tibble: 37 × 4
   f_occ          f_emp      f_class          n
   <fct>          <fct>      <fct>        <int>
 1 運輸           自営業     自営業          35
 2 運輸           正規雇用   マニュアル     629
 3 運輸           非正規雇用 非正規雇用      31
 4 運輸           <NA>       <NA>             2
 5 管理職         自営業     専門・管理職   106
 6 管理職         正規雇用   専門・管理職  1371
 7 管理職         非正規雇用 非正規雇用       6
 8 管理職         <NA>       <NA>            18
 9 技能・労務     自営業     自営業         411
10 技能・労務     正規雇用   マニュアル    2132
11 技能・労務     非正規雇用 非正規雇用     105
12 技能・労務     <NA>       <NA>            15
13 事務・営業     自営業     自営業          71
14 事務・営業     正規雇用   事務・販売職  3587
15 事務・営業     非正規雇用 非正規雇用      59
16 事務・営業     <NA>       <NA>             7
17 専門・技術     自営業     専門・管理職   397
18 専門・技術     正規雇用   専門・管理職  2860
19 専門・技術     非正規雇用 非正規雇用      42
20 専門・技術     <NA>       <NA>            25
21 農林漁業       自営業     自営業         102
22 農林漁業       正規雇用   マニュアル      42
23 農林漁業       非正規雇用 非正規雇用       1
24 農林漁業       <NA>       <NA>             2
25 販売・サービス 自営業     自営業         355
26 販売・サービス 正規雇用   事務・販売職  1144
27 販売・サービス 非正規雇用 非正規雇用      70
28 販売・サービス <NA>       <NA>             8
29 保安           自営業     自営業           1
30 保安           正規雇用   マニュアル     257
31 保安           非正規雇用 非正規雇用      12
32 保安           <NA>       <NA>             4
33 <NA>           自営業     自営業         243
34 <NA>           正規雇用   <NA>           805
35 <NA>           非正規雇用 非正規雇用      51
36 <NA>           無職       無職            98
37 <NA>           <NA>       <NA>         14742

# 母親の職業階層
d <- d |>
  mutate(
    m_class = case_when(
      # 専門・管理職
      m_occ == "管理職" & m_emp %in% c("正規雇用", "自営業") ~ "専門・管理職",
      m_occ == "専門・技術" & m_emp %in% c("正規雇用", "自営業") ~ "専門・管理職",
      # 事務・販売職（正規雇用のみ）
      m_occ == "事務・営業" & m_emp == "正規雇用" ~ "事務・販売職",
      m_occ == "販売・サービス" & m_emp == "正規雇用" ~ "事務・販売職",
      # 自営業（農林漁業の自営を含む）
      m_emp == "自営業" ~ "自営業",
      m_occ == "農林漁業" & m_emp == "自営業" ~ "自営業",
      # マニュアル（技能・運輸・保安で正規雇用、農林漁業の正規雇用を含む）
      m_occ %in% c("技能・労務", "運輸", "保安") & m_emp == "正規雇用" ~ "マニュアル",
      m_occ == "農林漁業" & m_emp == "正規雇用" ~ "マニュアル",
      # 非正規雇用（農林漁業の非正規を含む）
      m_occ == "農林漁業" & m_emp == "非正規雇用" ~ "非正規雇用",
      # 無職
      m_emp == "無職・専業主婦" ~ "無職",
      # 母親なし
      m_emp == "母親なし" ~ "母親なし",
      .default = NA
    ),
    m_class = factor(m_class, levels = c(
      "専門・管理職", "事務・販売職", "自営業",
      "マニュアル", "非正規雇用", "無職", "母親なし"
    ))
  )

# 分布を確認
d |> count(m_class)

# A tibble: 7 × 2
  m_class          n
  <fct>        <int>
1 専門・管理職  1281
2 事務・販売職  1462
3 自営業         715
4 マニュアル     122
5 非正規雇用      35
6 無職          4474
7 <NA>         21757

EGP分類について

EGP（Erikson-Goldthorpe-Portocarero）分類は国際比較に用いられる職業階層分類で、雇用関係（employment relationship）に基づく。本来は詳細な職業コードと雇用上の地位から分類するが、ここでは利用可能な変数から簡易的に作成している。厳密なEGP分類が必要な場合は、SSM調査などより詳細な職業情報を持つデータの利用を検討する。

coalesce()について

coalesce()は，複数の変数から最初の非欠損値を取得する．性別や親学歴のように時間で変化しない（または変化しにくい）変数は，どのWaveで回答されても同じ値のはずなので，この方法で欠損を補完できる．

# 各Waveの年度
# Wave1: 2015, Wave2: 2016, Wave3: 2017, Wave4: 2018,
# Wave5: 2019, Wave6: 2020, Wave7: 2021

# 在学中（学年1-12）の最初の観測から出生年度を計算
d <- d |> 
  mutate(
    # 在学中（1-12）のみを有効な学年として扱う
    w1学年_valid = if_else(w1学年 %in% 1:12, w1学年, NA),
    w2学年_valid = if_else(w2学年 %in% 1:12, w2学年, NA),
    w3学年_valid = if_else(w3学年 %in% 1:12, w3学年, NA),
    w4学年_valid = if_else(w4学年 %in% 1:12, w4学年, NA),
    w5学年_valid = if_else(w5学年 %in% 1:12, w5学年, NA),
    w6学年_valid = if_else(w6学年 %in% 1:12, w6学年, NA),
    w7学年_valid = if_else(w7学年 %in% 1:12, w7学年, NA),
    
    # 最初に観測された有効な学年
    first_grade = coalesce(w1学年_valid, w2学年_valid, w3学年_valid, 
                           w4学年_valid, w5学年_valid, w6学年_valid, 
                           w7学年_valid),
    # どのWaveで最初に有効な学年が観測されたか
    first_wave = case_when(
      !is.na(w1学年_valid) ~ 2015,
      !is.na(w2学年_valid) ~ 2016,
      !is.na(w3学年_valid) ~ 2017,
      !is.na(w4学年_valid) ~ 2018,
      !is.na(w5学年_valid) ~ 2019,
      !is.na(w6学年_valid) ~ 2020,
      !is.na(w7学年_valid) ~ 2021,
      .default = NA
    ),
    # 出生年度を計算
    birth_year = first_wave - first_grade - 6
  ) |> 
  # 作業用変数を削除
  select(-ends_with("_valid"))

d |> 
  count(first_wave)

# A tibble: 7 × 2
  first_wave     n
       <dbl> <int>
1       2015 18084
2       2016  1856
3       2017  1946
4       2018  1739
5       2019  1933
6       2020  2403
7       2021  1885

d |> 
  count(first_grade)

# A tibble: 12 × 2
   first_grade     n
   <dbl+lbl>   <int>
 1  1 [小１]   13111
 2  2 [小２]    1876
 3  3 [小３]    1719
 4  4 [小４]    1523
 5  5 [小５]    1416
 6  6 [小６]    1421
 7  7 [中１]    1488
 8  8 [中２]    1540
 9  9 [中３]    1511
10 10 [高１]    1437
11 11 [高２]    1420
12 12 [高３]    1384

d |> 
  count(birth_year)

# A tibble: 18 × 2
   birth_year     n
        <dbl> <int>
 1       1997  1382
 2       1998  1420
 3       1999  1421
 4       2000  1506
 5       2001  1541
 6       2002  1475
 7       2003  1420
 8       2004  1436
 9       2005  1465
10       2006  1628
11       2007  1591
12       2008  1837
13       2009  1857
14       2010  2016
15       2011  1990
16       2012  1870
17       2013  2106
18       2014  1885

# case_when：条件による分岐
d <- d |> 
  mutate(
    vocab_cat = case_when(
      vocab_w2_g3 >= 60 ~ "高",
      vocab_w2_g3 >= 40 ~ "中",
      vocab_w2_g3 < 40 ~ "低",
      .default = NA
    )
  )
d |> count(vocab_cat)

# A tibble: 4 × 2
  vocab_cat     n
  <chr>     <int>
1 中         2141
2 低          135
3 高         1357
4 <NA>      26213

.defaultの使い分け

• .default = NA：条件に合わないものは欠損値にしたいとき
• .default = 元の変数：一部だけ変更し，他は元のままにしたいとき
• .default = "その他"：条件に合わないものをまとめたいとき

2.4.7 時間とともに変化する変数

# ある年における学年（出生年度から計算）
d <- d |> 
  mutate(
    grade__2015 = 2015 - birth_year - 6,
    grade__2016 = 2016 - birth_year - 6,
    grade__2017 = 2017 - birth_year - 6,
    grade__2018 = 2018 - birth_year - 6,
    grade__2019 = 2019 - birth_year - 6,
    grade__2020 = 2020 - birth_year - 6,
    grade__2021 = 2021 - birth_year - 6
  )

# 確認
d |> count(birth_year, grade__2015, grade__2021)

# A tibble: 18 × 4
   birth_year grade__2015 grade__2021     n
        <dbl>       <dbl>       <dbl> <int>
 1       1997          12          18  1382
 2       1998          11          17  1420
 3       1999          10          16  1421
 4       2000           9          15  1506
 5       2001           8          14  1541
 6       2002           7          13  1475
 7       2003           6          12  1420
 8       2004           5          11  1436
 9       2005           4          10  1465
10       2006           3           9  1628
11       2007           2           8  1591
12       2008           1           7  1837
13       2009           0           6  1857
14       2010          -1           5  2016
15       2011          -2           4  1990
16       2012          -3           3  1870
17       2013          -4           2  2106
18       2014          -5           1  1885

# 学習塾へ
# 通塾変数の整理（年度ベース，0/1に変換）
d <- d |> 
  mutate(
    juku__2015 = case_match(ds0000000121_w1p,
                            1 ~ 1,  # 行っている
                            2 ~ 0,  # 行っていない
                            .default = NA),
    juku__2016 = case_match(ds0000000121_w2p,
                            1 ~ 1,
                            2 ~ 0,
                            .default = NA),
    juku__2017 = case_match(ds0000000121_w3p,
                            1 ~ 1,
                            2 ~ 0,
                            .default = NA),
    juku__2018 = case_match(ds0000000121_w4p,
                            1 ~ 1,
                            2 ~ 0,
                            .default = NA),
    juku__2019 = case_match(ds0000000121_w5p,
                            1 ~ 1,
                            2 ~ 0,
                            .default = NA),
    juku__2020 = case_match(ds0000000121_w6p,
                            1 ~ 1,
                            2 ~ 0,
                            .default = NA),
    juku__2021 = case_match(ds0000000121_w7p,
                            1 ~ 1,
                            2 ~ 0,
                            .default = NA)
  )

# 確認
d |> count(ds0000000121_w1p, juku__2015)

# A tibble: 4 × 3
  ds0000000121_w1p    juku__2015     n
  <dbl+lbl>                <dbl> <int>
1    1 [行っている]            1  5218
2    2 [行ってない]            0 11490
3 9999 [無回答・不明]         NA   102
4   NA                        NA 13036

# 年収の整理（年度ベースの変数名）
d <- d |> 
  mutate(
    income__2015 = if_else(ds0000000052_w1p %in% 1:10, 
                           as.numeric(ds0000000052_w1p), NA),
    income__2016 = if_else(ds0000000052_w2p %in% 1:10, 
                           as.numeric(ds0000000052_w2p), NA),
    income__2017 = if_else(ds0000000052_w3p %in% 1:10, 
                           as.numeric(ds0000000052_w3p), NA),
    income__2018 = if_else(ds0000000052_w4p %in% 1:10, 
                           as.numeric(ds0000000052_w4p), NA),
    income__2019 = if_else(ds0000000052_w5p %in% 1:10, 
                           as.numeric(ds0000000052_w5p), NA),
    income__2020 = if_else(ds0000000052_w6p %in% 1:10, 
                           as.numeric(ds0000000052_w6p), NA),
    income__2021 = if_else(ds0000000052_w7p %in% 1:10, 
                           as.numeric(ds0000000052_w7p), NA)
  )

年度ベースの変数名の利点

birth_yearと組み合わせて年齢が計算しやすい（年度 - birth_year） pivot_longer()で変換後，year変数がそのまま年度になる 直感的に「いつのデータか」が分かりやすい アンダースコア2つ（）で区切ることで，names_sep = “”での分割が容易

d <- d |> 
  mutate(across(income__2015:income__2021,
                ~ case_match(.x,
      1 ~ 100,
      2 ~ 250,
      3 ~ 350,
      4 ~ 450,
      5 ~ 550,
      6 ~ 700,
      7 ~ 900,
      8 ~ 1250,
      9 ~ 1750,
      10 ~ 2500,  # 上限なしのため仮の値
      .default = NA
    )
  ))

d |> count(income__2015)

# A tibble: 11 × 2
   income__2015     n
          <dbl> <int>
 1          100   414
 2          250   756
 3          350  1444
 4          450  1973
 5          550  2545
 6          700  3597
 7          900  2242
 8         1250  1627
 9         1750   242
10         2500   119
11           NA 14887

d |> count(income__2021)

# A tibble: 11 × 2
   income__2021     n
          <dbl> <int>
 1          100   347
 2          250   594
 3          350   912
 4          450  1387
 5          550  2011
 6          700  3535
 7          900  2375
 8         1250  2024
 9         1750   389
10         2500   160
11           NA 16112

# 将来の希望進学段階（子ども）→ 教育年数に変換
# 値ラベル: 1=中学まで(9年), 2=高校まで(12年), 3=専門学校まで(14年),
#          4=短大まで(14年), 5=大学まで(16年), 6=大学院まで(18年)
# 7=その他, 8=まだ決めていない は欠損扱い

d <- d |>
  mutate(
    # 教育年数に変換
    edu_asp__2015 = case_match(as.numeric(ds0000000198_w1c),
                               1 ~ 9, 2 ~ 12, 3 ~ 14, 4 ~ 14, 5 ~ 16, 6 ~ 18,
                               .default = NA),
    edu_asp__2016 = case_match(as.numeric(ds0000000198_w2c),
                               1 ~ 9, 2 ~ 12, 3 ~ 14, 4 ~ 14, 5 ~ 16, 6 ~ 18,
                               .default = NA),
    edu_asp__2017 = case_match(as.numeric(ds0000000198_w3c),
                               1 ~ 9, 2 ~ 12, 3 ~ 14, 4 ~ 14, 5 ~ 16, 6 ~ 18,
                               .default = NA),
    edu_asp__2018 = case_match(as.numeric(ds0000000198_w4c),
                               1 ~ 9, 2 ~ 12, 3 ~ 14, 4 ~ 14, 5 ~ 16, 6 ~ 18,
                               .default = NA),
    edu_asp__2019 = case_match(as.numeric(ds0000000198_w5c),
                               1 ~ 9, 2 ~ 12, 3 ~ 14, 4 ~ 14, 5 ~ 16, 6 ~ 18,
                               .default = NA),
    edu_asp__2020 = case_match(as.numeric(ds0000000198_w6c),
                               1 ~ 9, 2 ~ 12, 3 ~ 14, 4 ~ 14, 5 ~ 16, 6 ~ 18,
                               .default = NA),
    edu_asp__2021 = case_match(as.numeric(ds0000000198_w7c),
                               1 ~ 9, 2 ~ 12, 3 ~ 14, 4 ~ 14, 5 ~ 16, 6 ~ 18,
                               .default = NA),

    # 大学進学希望ダミー（大学・大学院=1, それ以外=0）
    univ_asp__2015 = case_when(ds0000000198_w1c %in% 5:6 ~ 1,
                               ds0000000198_w1c %in% 1:4 ~ 0,
                               .default = NA),
    univ_asp__2016 = case_when(ds0000000198_w2c %in% 5:6 ~ 1,
                               ds0000000198_w2c %in% 1:4 ~ 0,
                               .default = NA),
    univ_asp__2017 = case_when(ds0000000198_w3c %in% 5:6 ~ 1,
                               ds0000000198_w3c %in% 1:4 ~ 0,
                               .default = NA),
    univ_asp__2018 = case_when(ds0000000198_w4c %in% 5:6 ~ 1,
                               ds0000000198_w4c %in% 1:4 ~ 0,
                               .default = NA),
    univ_asp__2019 = case_when(ds0000000198_w5c %in% 5:6 ~ 1,
                               ds0000000198_w5c %in% 1:4 ~ 0,
                               .default = NA),
    univ_asp__2020 = case_when(ds0000000198_w6c %in% 5:6 ~ 1,
                               ds0000000198_w6c %in% 1:4 ~ 0,
                               .default = NA),
    univ_asp__2021 = case_when(ds0000000198_w7c %in% 5:6 ~ 1,
                               ds0000000198_w7c %in% 1:4 ~ 0,
                               .default = NA)
  )

# 確認
d |> count(ds0000000198_w1c, edu_asp__2015, univ_asp__2015)

# A tibble: 10 × 4
   ds0000000198_w1c                  edu_asp__2015 univ_asp__2015     n
   <dbl+lbl>                                 <dbl>          <dbl> <int>
 1    1 [中学校まで]                             9              0    30
 2    2 [高校まで]                              12              0  1037
 3    3 [専門学校・各種学校まで]                14              0  1022
 4    4 [短期大学まで]                          14              0   353
 5    5 [大学（四年制、六年制）まで]            16              1  6448
 6    6 [大学院まで]                            18              1   772
 7    7 [その他]                                NA             NA    34
 8    8 [まだ決めていない]                      NA             NA  2187
 9 9999 [無回答・不明]                          NA             NA   119
10   NA                                         NA             NA 17844

将来の希望進学段階の変数

• edu_asp__年度: 希望進学段階を教育年数に変換（中学9年〜大学院18年）
• univ_asp__年度: 大学進学希望ダミー（大学・大学院希望=1, その他=0）

元の値	進学段階	教育年数
1	中学まで	9年
2	高校まで	12年
3	専門学校まで	14年
4	短大まで	14年
5	大学まで	16年
6	大学院まで	18年

「その他」「まだ決めていない」は分析上の解釈が難しいため欠損として処理している．

2.4.8 summarise()：集計

# 全体の平均と標準偏差
d |> 
  summarise(
    mean_vocab = mean(vocab_w2_g3, na.rm = TRUE),
    sd_vocab = sd(vocab_w2_g3, na.rm = TRUE),
    n = sum(!is.na(vocab_w2_g3))
  )

# A tibble: 1 × 3
  mean_vocab sd_vocab     n
       <dbl>    <dbl> <int>
1       57.4     11.4  3633

# グループ別の集計（.by引数を使用）
d |> 
  summarise(
    mean_vocab = mean(vocab_w2_g3, na.rm = TRUE),
    sd_vocab = sd(vocab_w2_g3, na.rm = TRUE),
    n = sum(!is.na(vocab_w2_g3)),
    .by = sex
  )

# A tibble: 3 × 4
  sex   mean_vocab sd_vocab     n
  <chr>      <dbl>    <dbl> <int>
1 男子        57.6     11.6  1727
2 女子        57.2     11.3  1877
3 <NA>        56.5     11.3    29

欠損値への対応

mean()やsd()は，データに欠損値（NA）が含まれると結果がNAになる．na.rm = TRUEを指定することで欠損値を除いて計算する．また，n()は欠損値を含む全ケース数を返すため，実際に計算に使用されたケース数を知りたい場合はsum(!is.na(変数名))を使用する．

2.5 自作関数の作成

繰り返し使う処理は，自作関数にまとめると便利である．

2.5.1 関数の基本構造

関数名 <- function(引数1, 引数2, ...) {
  # 処理内容
  return(結果)  # returnは省略可（最後の式が返される）
}

2.5.2 95%信頼区間付きの記述統計関数

記述統計と信頼区間を一度に計算する関数を作成する．

# 95%信頼区間付きの記述統計を計算する関数
# データフレームと変数名を受け取る（tidyverse風）
desc_with_ci <- function(data, var) {
  data |>
    dplyr::summarise(
      Mean = mean({{ var }}, na.rm = TRUE),
      SD = sd({{ var }}, na.rm = TRUE),
      n = sum(!is.na({{ var }})),
      SE = SD / sqrt(n),
      CI_lower = Mean + qt(0.025, df = n - 1) * SE,
      CI_upper = Mean + qt(0.975, df = n - 1) * SE
    )
}

# 使用例
desc_with_ci(d, vocab_w2_g3)

# A tibble: 1 × 6
   Mean    SD     n    SE CI_lower CI_upper
  <dbl> <dbl> <int> <dbl>    <dbl>    <dbl>
1  57.4  11.4  3633 0.190     57.0     57.8

qt()による信頼区間の計算

上記のコードで使われているqt(0.025, df = n - 1)とqt(0.975, df = n - 1)は、t分布の分位関数である。自由度n - 1のt分布において、累積確率が0.025と0.975になる点（両側5%の臨界値）を返す。サンプルサイズが大きい場合、これらの値は正規分布の±1.96に近づく。

{ }（curly-curly）について

{ var }は，関数の引数として渡された変数名をそのままdplyr関数内で使えるようにする記法である．これによりdesc_with_ci(d, vocab_w2_g3)のように，クォートなしで変数名を指定できる．

2.5.3 グループ別に適用する

.by引数を追加すれば，グループ別の集計も簡単にできる．

# グループ別に対応した関数
desc_with_ci_by <- function(data, var, by = NULL) {
  data |>
    dplyr::summarise(
      Mean = mean({{ var }}, na.rm = TRUE),
      SD = sd({{ var }}, na.rm = TRUE),
      n = sum(!is.na({{ var }})),
      SE = SD / sqrt(n),
      CI_lower = Mean + qt(0.025, df = n - 1) * SE,
      CI_upper = Mean + qt(0.975, df = n - 1) * SE,
      .by = {{ by }}
    )
}

# 父親の学歴別に適用
d |>
  drop_na(f_edu) |>
  desc_with_ci_by(vocab_w2_g3, by = f_edu)

# A tibble: 3 × 7
  f_edu       Mean    SD     n    SE CI_lower CI_upper
  <fct>      <dbl> <dbl> <int> <dbl>    <dbl>    <dbl>
1 専門・短大  56.1 10.6    454 0.496     55.1     57.0
2 中学・高校  56.1  9.97   899 0.333     55.5     56.8
3 大学・院    58.4 12.0   1650 0.295     57.9     59.0

自作関数の利点

• 同じ処理を何度も書かなくてよい
• コードが読みやすくなる
• 修正が一箇所で済む（DRY原則：Don’t Repeat Yourself）

2.6 記述統計

2.6.1 量的変数の要約

# 語彙力スコアの記述統計
d |> 
  summarise(
    Mean = mean(vocab_w2_g3, na.rm = TRUE),
    SD = sd(vocab_w2_g3, na.rm = TRUE),
    Min = min(vocab_w2_g3, na.rm = TRUE),
    Max = max(vocab_w2_g3, na.rm = TRUE),
    Median = median(vocab_w2_g3, na.rm = TRUE),
    n = sum(!is.na(vocab_w2_g3))
  )

# A tibble: 1 × 6
   Mean    SD   Min   Max Median     n
  <dbl> <dbl> <dbl> <dbl>  <dbl> <int>
1  57.4  11.4  27.4  105.   56.5  3633

2.6.2 グループ別の記述統計

# 性別ごとの記述統計
d |> 
  summarise(
    Mean = mean(vocab_w2_g3, na.rm = TRUE),
    SD = sd(vocab_w2_g3, na.rm = TRUE),
    n = sum(!is.na(vocab_w2_g3)),
    .by = sex
  )

# A tibble: 3 × 4
  sex    Mean    SD     n
  <chr> <dbl> <dbl> <int>
1 男子   57.6  11.6  1727
2 女子   57.2  11.3  1877
3 <NA>   56.5  11.3    29

# 親の学歴別の記述統計
d |> 
  summarise(
    Mean = mean(vocab_w2_g3, na.rm = TRUE),
    SD = sd(vocab_w2_g3, na.rm = TRUE),
    n = sum(!is.na(vocab_w2_g3)),
    .by = f_edu
  )

# A tibble: 4 × 4
  f_edu       Mean    SD     n
  <fct>      <dbl> <dbl> <int>
1 専門・短大  56.1 10.6    454
2 <NA>        57.5 12.3    630
3 中学・高校  56.1  9.97   899
4 大学・院    58.4 12.0   1650

# 親の学歴別の記述統計
d |> 
  summarise(
    Mean = mean(income__2015, na.rm = TRUE),
    SD = sd(income__2015, na.rm = TRUE),
    n = sum(!is.na(income__2015)),
    .by = f_edu
  ) |> 
  arrange(f_edu)

# A tibble: 4 × 4
  f_edu       Mean    SD     n
  <fct>      <dbl> <dbl> <int>
1 中学・高校  585.  257.  4217
2 専門・短大  599.  266.  2147
3 大学・院    829.  380.  7067
4 <NA>        462.  326.  1528

# 親の学歴別の記述統計
d |> 
  summarise(
    Mean = mean(income__2015, na.rm = TRUE),
    SD = sd(income__2015, na.rm = TRUE),
    n = sum(!is.na(income__2015)),
    .by = c(f_edu, m_edu)
  ) |> 
  arrange(f_edu, m_edu)

# A tibble: 16 × 5
   f_edu      m_edu       Mean    SD     n
   <fct>      <fct>      <dbl> <dbl> <int>
 1 中学・高校 中学・高校  557.  247.  1927
 2 中学・高校 専門・短大  600.  255.  1774
 3 中学・高校 大学・院    650.  289.   446
 4 中学・高校 <NA>        542.  245.    70
 5 専門・短大 中学・高校  548.  252.   565
 6 専門・短大 専門・短大  607.  258.  1213
 7 専門・短大 大学・院    666.  306.   327
 8 専門・短大 <NA>        539.  195.    42
 9 大学・院   中学・高校  704.  331.   901
10 大学・院   専門・短大  795.  353.  3005
11 大学・院   大学・院    902.  406.  3047
12 大学・院   <NA>        745.  297.   114
13 <NA>       中学・高校  288.  241.   351
14 <NA>       専門・短大  421.  302.   527
15 <NA>       大学・院    520.  381.   234
16 <NA>       <NA>        628.  295.   416

2.7 度数分布表

janitorパッケージのtabyl()が便利である．

# 性別の度数分布
d |> tabyl(sex)

  sex     n    percent valid_percent
 女子 14607 0.48941232     0.4986856
 男子 14684 0.49199223     0.5013144
 <NA>   555 0.01859546            NA

# 父親の学歴の度数分布
d |> tabyl(f_edu)

      f_edu     n   percent valid_percent
 中学・高校  7326 0.2454600     0.2958088
 専門・短大  4105 0.1375394     0.1657514
   大学・院 13335 0.4467935     0.5384398
       <NA>  5080 0.1702071            NA

2.8 相関係数

2つの連続変数間の線形な関連の強さを表す指標である．-1から1の値をとり，絶対値が大きいほど関連が強い．

# 2変数間の相関係数
cor(d$vocab_w2_g3, d$vocab_w2_g9, use = "complete.obs")

[1] 0.5120338

use引数について

cor()関数のuse引数は欠損値の扱いを指定する．

• "everything": 欠損があるとNAを返す（デフォルト）
• "complete.obs": 両方の変数で欠損のないケースのみ使用
• "pairwise.complete.obs": 変数のペアごとに欠損のないケースを使用

2.8.1 相関係数の検定

相関係数が統計的に有意かどうかを検定する．

# 相関係数の検定
cor.test(d$vocab_w2_g3, d$vocab_w2_g9)

    Pearson's product-moment correlation

data:  d$vocab_w2_g3 and d$vocab_w2_g9
t = 19.417, df = 1061, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.4662615 0.5550719
sample estimates:
      cor 
0.5120338 

2.8.2 相関行列

複数の変数間の相関を一度に確認したい場合は，相関行列を作成する．

# 分析に使う変数を選択
d_cor <- d |>
  select(vocab_w2_g3, vocab_w2_g9, edu_asp__2016, income__2016) |>
  drop_na()

# 相関行列
cor(d_cor)

              vocab_w2_g3 vocab_w2_g9 edu_asp__2016 income__2016
vocab_w2_g3     1.0000000  0.48932554     0.1908636   0.11962830
vocab_w2_g9     0.4893255  1.00000000     0.1901056   0.05715703
edu_asp__2016   0.1908636  0.19010563     1.0000000   0.20493341
income__2016    0.1196283  0.05715703     0.2049334   1.00000000

# 見やすく丸める
cor(d_cor) |> round(2)

              vocab_w2_g3 vocab_w2_g9 edu_asp__2016 income__2016
vocab_w2_g3          1.00        0.49          0.19         0.12
vocab_w2_g9          0.49        1.00          0.19         0.06
edu_asp__2016        0.19        0.19          1.00         0.20
income__2016         0.12        0.06          0.20         1.00

2.8.3 相関行列の可視化

corrplotパッケージを使うと相関行列を視覚的に確認できる．

# パッケージのインストール（初回のみ）
install.packages("corrplot")

library(corrplot)
corrplot(cor(d_cor), method = "color", type = "upper",
         addCoef.col = "black", number.cex = 0.8,
         tl.col = "black", tl.srt = 45)

相関係数の解釈

相関係数はあくまで2変数間の線形な関連を表す．

• 因果関係を意味しない
• 非線形な関連は捉えられない
• 第3の変数（交絡因子）の影響を受ける

社会科学では，相関係数の絶対値が0.3程度でも「中程度の関連」と解釈されることが多い．

2.9 ggplotによる可視化

ggplot2パッケージを使って図を作成する．

信頼区間を示すことの重要性

可視化において信頼区間（エラーバー）を示すことは非常に重要である．平均値や回帰係数などの点推定値だけを示すと，推定の不確実性が伝わらない．信頼区間を示すことで：

• 推定値の精度（サンプルサイズや分散の影響）が視覚的に分かる
• グループ間の差が統計的に有意かどうかの目安になる
• 読者が結果の解釈を適切に行える

本セクションでは，geom_errorbar()やgeom_smooth(se = TRUE)を使って信頼区間を表示する方法を紹介する．

注意：2つのグループの95%信頼区間が重なっていても，その差が統計的に有意でないとは限らない．信頼区間の重複は目安に過ぎず，正確な検定にはt検定や回帰分析を行う必要がある．

Macでの日本語表示

Macではggplot2の日本語が文字化けすることがある．以下のいずれかで対処できる．

方法1：テーマで指定

theme_set(theme_minimal(base_family = "HiraKakuProN-W3"))

方法2：個別の図で指定

ggplot(...) +
  ... +
  theme_minimal(base_family = "HiraKakuProN-W3")

方法3：raggパッケージを使用

install.packages("ragg")

RStudioの設定で「Tools」→「Global Options」→「General」→「Graphics」→「Backend」を「AGG」に変更する．

ここでは次のようなテーマを設定する．

# MACの場合
theme_set(theme_minimal(base_family = "HiraKakuProN-W3"))
# Windowsの場合 #をはずす
# theme_set(theme_minimal())

2.9.1 ヒストグラム

d |> 
  ggplot(aes(x = vocab_w2_g3)) +
  geom_histogram(binwidth = 5, 
                 fill = "steelblue",
                 color = "white") +
  labs(x = "語彙力スコア（小3基準偏差値）", 
       y = "度数", 
       title = "語彙力スコアの分布")

2.9.2 箱ひげ図

d |> 
  drop_na(sex) |> 
  ggplot(aes(x = sex, y = vocab_w2_g3)) +
  geom_boxplot(fill = "lightblue") +
  labs(x = "性別", y = "語彙力スコア")

2.9.3 棒グラフ

棒グラフには2つの方法がある．

方法1：geom_bar()（集計前のデータ）

geom_bar()は自動的に度数を集計する．

d |> 
  drop_na(f_edu) |> 
  ggplot(aes(x = f_edu)) +
  geom_bar(fill = "steelblue") +
  labs(x = "父親の学歴", y = "度数")

方法2：geom_col()（集計後のデータ）

事前に集計したデータを使う場合はgeom_col()を使用する．

# 集計
d_count <- d |> 
  drop_na(f_edu) |> 
  count(f_edu)
d_count

# A tibble: 3 × 2
  f_edu          n
  <fct>      <int>
1 中学・高校  7326
2 専門・短大  4105
3 大学・院   13335

# 棒グラフ
d_count |> 
  ggplot(aes(x = f_edu, y = n)) +
  geom_col(fill = "steelblue") +
  labs(x = "父親の学歴", y = "度数")

平均値と信頼区間の棒グラフ

集計後のデータを使えば，信頼区間を追加できる．

# 父親の学歴別の平均と信頼区間を集計
d_summary <- d |> 
  drop_na(f_edu) |> 
  summarise(
    Mean = mean(vocab_w2_g3, na.rm = TRUE),
    SD = sd(vocab_w2_g3, na.rm = TRUE),
    n = sum(!is.na(vocab_w2_g3)),
    SE = SD / sqrt(n),
    ll = Mean + qt(0.025, df = n - 1) * SE,
    ul = Mean + qt(0.975, df = n - 1) * SE,
    .by = f_edu
  )
d_summary

# A tibble: 3 × 7
  f_edu       Mean    SD     n    SE    ll    ul
  <fct>      <dbl> <dbl> <int> <dbl> <dbl> <dbl>
1 専門・短大  56.1 10.6    454 0.496  55.1  57.0
2 中学・高校  56.1  9.97   899 0.333  55.5  56.8
3 大学・院    58.4 12.0   1650 0.295  57.9  59.0

# 棒グラフ＋エラーバー
d_summary |> 
  ggplot(aes(x = f_edu, y = Mean, ymin = ll, ymax = ul)) +
  geom_col(fill = "steelblue") +
  geom_errorbar(width = 0.2) +
  labs(x = "父親の学歴", 
       y = "語彙力スコア（平均値と95%信頼区間）") 

2.9.4 散布図

d |> 
  ggplot(aes(x = vocab_w2_g3, 
             y = vocab_w2_g9)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = TRUE) +
  labs(x = "語彙力スコア", y = "読解力スコア") 

2.9.5 グループ別の図

d |>
  drop_na(sex) |>
  ggplot(aes(x = vocab_w2_g3,
             y = vocab_w2_g9,
             color = sex, shape = sex)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = FALSE, linetype = "solid") +
  labs(x = "語彙力スコア",
       y = "読解力スコア", color = "性別", shape = "性別") +
  scale_color_viridis_d(option = "C", end = .8)

2.9.6 ファセット（パネル分割）

d |> 
  drop_na(sex) |> 
  ggplot(aes(x = vocab_w2_g3)) +
  geom_histogram(binwidth = 5, 
                 fill = "steelblue", 
                 color = "white") +
  facet_wrap(~ sex) +
  labs(x = "語彙力スコア", y = "度数") 

2.10 平均値と信頼区間の図示

# 父親の学歴別の平均値と信頼区間
tab_vocab_by_fedu <- d |> 
  drop_na(f_edu) |> 
  summarise(
    Mean = mean(vocab_w2_g3, na.rm = TRUE),
    SD = sd(vocab_w2_g3, na.rm = TRUE),
    n = sum(!is.na(vocab_w2_g3)),
    SE = SD / sqrt(n),
    ll = Mean + qt(0.025, df = n - 1) * SE,
    ul = Mean + qt(0.975, df = n - 1) * SE,
    .by = f_edu
  )
tab_vocab_by_fedu

# A tibble: 3 × 7
  f_edu       Mean    SD     n    SE    ll    ul
  <fct>      <dbl> <dbl> <int> <dbl> <dbl> <dbl>
1 専門・短大  56.1 10.6    454 0.496  55.1  57.0
2 中学・高校  56.1  9.97   899 0.333  55.5  56.8
3 大学・院    58.4 12.0   1650 0.295  57.9  59.0

# 図示
tab_vocab_by_fedu |> 
  ggplot(aes(x = f_edu, y = Mean, ymin = ll, ymax = ul)) +
  geom_pointrange() +
  labs(x = "父親の学歴", 
       y = "語彙力スコア（平均値と95%信頼区間）")

3 実践編

3.1 クロス表分析

3.1.1 クロス表の作成

# 性別×親学歴のクロス表
d |> 
  drop_na(sex, f_edu) |> 
  tabyl(sex, f_edu)

  sex 中学・高校 専門・短大 大学・院
 女子       3626       2039     6658
 男子       3688       2064     6648

# 行パーセント
d |> 
  drop_na(sex, f_edu) |> 
  tabyl(sex, f_edu) |> 
  adorn_percentages("row") |> 
  adorn_pct_formatting(digits = 1)

  sex 中学・高校 専門・短大 大学・院
 女子      29.4%      16.5%    54.0%
 男子      29.7%      16.6%    53.6%

# 列パーセント
d |> 
  drop_na(sex, f_edu) |> 
  tabyl(sex, f_edu) |> 
  adorn_percentages("col") |> 
  adorn_pct_formatting(digits = 1)

  sex 中学・高校 専門・短大 大学・院
 女子      49.6%      49.7%    50.0%
 男子      50.4%      50.3%    50.0%

3.1.2 カイ二乗検定

# クロス表を作成（方法1：table関数）
d_complete <- d |>
  drop_na(sex, f_edu)
tab <- table(d_complete$sex, d_complete$f_edu)
tab

      
       中学・高校 専門・短大 大学・院
  女子       3626       2039     6658
  男子       3688       2064     6648

# クロス表を作成（方法2：xtabs関数）
# フォーミュラで変数を指定できる
tab <- xtabs(~ sex + f_edu, data = d_complete)
tab

      f_edu
sex    中学・高校 専門・短大 大学・院
  女子       3626       2039     6658
  男子       3688       2064     6648

# カイ二乗検定
chisq.test(tab)

    Pearson's Chi-squared test

data:  tab
X-squared = 0.4456, df = 2, p-value = 0.8003

3.2 t検定

2群の平均値の差を検定するにはt.test()関数を使用する．

# 性別による読解力スコアの差（2群比較）
t.test(vocab_w2_g9 ~ sex, data = d)

    Welch Two Sample t-test

data:  vocab_w2_g9 by sex
t = -0.30489, df = 1106.3, p-value = 0.7605
alternative hypothesis: true difference in means between group 女子 and group 男子 is not equal to 0
95 percent confidence interval:
 -1.453686  1.062668
sample estimates:
mean in group 女子 mean in group 男子 
          51.62289           51.81840 

# 等分散を仮定するt検定（Student's t-test）
t.test(vocab_w2_g9 ~ sex, data = d, var.equal = TRUE)

    Two Sample t-test

data:  vocab_w2_g9 by sex
t = -0.30504, df = 1113, p-value = 0.7604
alternative hypothesis: true difference in means between group 女子 and group 男子 is not equal to 0
95 percent confidence interval:
 -1.453076  1.062058
sample estimates:
mean in group 女子 mean in group 男子 
          51.62289           51.81840 

t検定の種類

• Welchのt検定（デフォルト）：等分散を仮定しない．ほとんどの場合こちらを使用
• Studentのt検定（var.equal = TRUE）：等分散を仮定．2群の分散が等しいと確信できる場合のみ
• 対応のあるt検定（paired = TRUE）：同一個人の前後比較など

Rの確率分布関数

Rでは各確率分布に対して4種類の関数が用意されている。t分布（t）、正規分布（norm）、カイ二乗分布（chisq）、F分布（f）などに共通の命名規則である。

接頭辞	関数名例	機能	用途
d	dt(), dnorm()	密度関数（確率密度）	分布の形状を描画
p	pt(), pnorm()	分布関数（累積確率）	p値の計算
q	qt(), qnorm()	分位関数（パーセンタイル点）	臨界値・信頼区間の計算
r	rt(), rnorm()	乱数生成	シミュレーション

# 例：t分布（自由度10）
dt(2, df = 10)       # t=2での確率密度
pt(2, df = 10)       # t≤2の累積確率（片側p値の補数）
qt(0.975, df = 10)   # 累積確率0.975に対応するt値（両側5%の臨界値）
rt(100, df = 10)     # t分布に従う乱数を100個生成

# 例：正規分布
qnorm(0.975)         # 1.96（95%信頼区間の係数）
pnorm(1.96)          # 0.975

3.2.1 順位変換とノンパラメトリック検定

データに外れ値が多い場合や正規性の仮定が満たされない場合，順位（rank）に変換してから分析する方法がある．

# 順位への変換
d <- d |>
  mutate(vocab_rank = rank(vocab_w2_g9, na.last = "keep"))

# 順位変換後のt検定（実質的にWilcoxon検定と同等）
t.test(vocab_rank ~ sex, data = d)

    Welch Two Sample t-test

data:  vocab_rank by sex
t = 0.095429, df = 1107.4, p-value = 0.924
alternative hypothesis: true difference in means between group 女子 and group 男子 is not equal to 0
95 percent confidence interval:
 -36.07326  39.76155
sample estimates:
mean in group 女子 mean in group 男子 
          559.0617           557.2176 

# Wilcoxon順位和検定（Mann-Whitney U検定）
# 順位に基づくノンパラメトリック検定
wilcox.test(vocab_w2_g9 ~ sex, data = d)

    Wilcoxon rank sum test with continuity correction

data:  vocab_w2_g9 by sex
W = 155770, p-value = 0.9229
alternative hypothesis: true location shift is not equal to 0

順位変換の利点

• 外れ値の影響を軽減できる
• 正規性の仮定が不要
• 元のスケールに依存しない比較が可能
• rank()で順位変換し，通常のt検定・回帰分析を適用する方法も実用的

Rank-Rank Correlation（順位-順位相関）

所得の世代間移動研究では，親の所得順位と子の所得順位の相関（rank-rank correlation）がよく用いられる（Chetty et al., 2014など）．所得を順位に変換することで：

• 所得分布の非正規性や外れ値の影響を軽減
• 異なる時代・地域間での比較が容易
• 「親が所得分布の下位25%にいる子が，上位25%に到達する確率」のような解釈が可能

順位変換の関数：

• rank(): 順位を返す（1, 2, 3, …）。ties.method引数でタイ（同順位）の扱いを指定（デフォルトは"average"で平均順位）
• percent_rank(): 0〜1のパーセンタイル順位を返す（(rank - 1) / (n - 1)）。タイは平均順位で処理
• dplyr::ntile(): n分位（四分位ならntile(x, 4)）

欠損値の扱い：rank()はna.last引数で制御（デフォルトはTRUEで末尾に配置）。percent_rank()は欠損値に対してNAを返す。

（0.5, 0.5）を通る性質：パーセンタイル順位の平均は常に0.5になるため，rank-rank回帰は必ず（0.5, 0.5）を通る。つまり「50パーセンタイルの親の子どもは，平均的に50パーセンタイルになる」という解釈が可能。

# 親子の所得順位相関の例
d |>
  mutate(
    parent_income_rank = percent_rank(parent_income),
    child_income_rank = percent_rank(child_income)
  ) |>
  summarise(rank_rank_cor = cor(parent_income_rank, child_income_rank,
                                 use = "complete.obs"))

3.3 分散分析（ANOVA）

3群以上の平均値の差を検定するには分散分析を使用する．

# 父親学歴（3群）による読解力スコアの差
model_aov <- aov(vocab_w2_g9 ~ f_edu, data = d)
summary(model_aov)

             Df Sum Sq Mean Sq F value   Pr(>F)    
f_edu         2   4233  2116.5   19.12 7.31e-09 ***
Residuals   914 101156   110.7                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
28929 個の観測値が欠損のため削除されました

# ANOVAで有意差があった場合，どの群間に差があるかを調べる（Tukeyの多重比較）
TukeyHSD(model_aov)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = vocab_w2_g9 ~ f_edu, data = d)

$f_edu
                           diff       lwr      upr     p adj
専門・短大-中学・高校 -1.005460 -3.849226 1.838306 0.6844689
大学・院-中学・高校    4.060963  2.181938 5.939988 0.0000014
大学・院-専門・短大    5.066423  2.470574 7.662271 0.0000156

3.3.1 効果量

p値は統計的有意性を示すが、効果の大きさは示さない。effectsizeパッケージで効果量を計算できる。

library(effectsize)

# η²（イータ二乗）：説明された分散の割合
eta_squared(model_aov)

# Effect Size for ANOVA

Parameter | Eta2 |       95% CI
-------------------------------
f_edu     | 0.04 | [0.02, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

# ω²（オメガ二乗）：より保守的な推定値（母集団への一般化を考慮）
omega_squared(model_aov)

# Effect Size for ANOVA

Parameter | Omega2 |       95% CI
---------------------------------
f_edu     |   0.04 | [0.02, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

効果量の解釈の目安（Cohen, 1988）

η²やω²の解釈の目安：

• 小: 0.01
• 中: 0.06
• 大: 0.14

ただし、分野や文脈によって「大きな効果」の基準は異なる。教育研究では0.01程度でも実質的に重要な場合がある。

多重比較の問題

複数の群を比較する場合，検定の回数が増えるほど偶然有意になる確率（第1種の過誤）が高くなる．3群の場合，すべての組み合わせ（3回）を個別にt検定すると，少なくとも1つが偶然有意になる確率は約14%に上昇する．TukeyのHSD法などの多重比較法は，この問題を補正する．

ANOVAと回帰分析の関係

カテゴリ変数を説明変数とする回帰分析は，ANOVAと本質的に同じ分析である．回帰分析では参照カテゴリとの差として係数が出力され，ANOVAではF検定で全体の差を検定する．

3.4 回帰分析

回帰分析は，説明変数\(X\)を条件としたときの従属変数\(Y\)の期待値\(E[Y|X]\)を線形予測子でモデル化する手法である．例えば，語彙力スコアを条件としたときの読解力スコアの期待値\(E[\text{読解力}|\text{語彙力}]\)をモデル化する．Rではlm()関数で線形回帰モデルを推定できる．

回帰係数の解釈について

回帰分析で得られる係数は，他の変数を統制した上での関連（条件付き関連）を示すものである．観察データから得られた係数を直ちに因果効果として解釈することには注意が必要である．因果効果として解釈するためには，交絡変数の統制，選択バイアスの考慮，時間的順序の確認など，追加の検討が必要となる．

3.4.1 単回帰分析

1つの説明変数で従属変数を予測するモデルである．

# 語彙力と読解力の関連
model1 <- lm(vocab_w2_g9 ~ vocab_w2_g3,
             data = d)
summary(model1)

Call:
lm(formula = vocab_w2_g9 ~ vocab_w2_g3, data = d)

Residuals:
    Min      1Q  Median      3Q     Max 
-30.545  -5.770  -0.146   4.877  37.352 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 18.32143    1.75282   10.45   <2e-16 ***
vocab_w2_g3  0.52007    0.02678   19.42   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.215 on 1061 degrees of freedom
  (28783 個の観測値が欠損のため削除されました)
Multiple R-squared:  0.2622,    Adjusted R-squared:  0.2615 
F-statistic:   377 on 1 and 1061 DF,  p-value: < 2.2e-16

3.4.2 結果の整理：broomパッケージ

broomパッケージには，モデルの結果をtibble形式で取得する3つの関数がある．

関数	内容	行数
tidy()	係数・標準誤差・p値など	説明変数の数
glance()	R²・AICなどモデル全体の指標	1行
augment()	予測値・残差など観測ごとの値	観測数

# tidy(): 係数の情報
tidy(model1)

# A tibble: 2 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)   18.3      1.75        10.5 2.08e-24
2 vocab_w2_g3    0.520    0.0268      19.4 4.23e-72

# tidy()に信頼区間を追加
tidy(model1, conf.int = TRUE)

# A tibble: 2 × 7
  term        estimate std.error statistic  p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
1 (Intercept)   18.3      1.75        10.5 2.08e-24   14.9      21.8  
2 vocab_w2_g3    0.520    0.0268      19.4 4.23e-72    0.468     0.573

# glance(): モデル全体の適合度指標
glance(model1)

# A tibble: 1 × 12
  r.squared adj.r.squared sigma statistic  p.value    df logLik   AIC   BIC
      <dbl>         <dbl> <dbl>     <dbl>    <dbl> <dbl>  <dbl> <dbl> <dbl>
1     0.262         0.261  9.21      377. 4.23e-72     1 -3868. 7742. 7757.
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

# augment(): 観測ごとの予測値・残差
augment(model1) |> head()

# A tibble: 6 × 9
  .rownames vocab_w2_g9 vocab_w2_g3 .fitted .resid    .hat .sigma   .cooksd
  <chr>           <dbl>       <dbl>   <dbl>  <dbl>   <dbl>  <dbl>     <dbl>
1 15784            70.5        71.4    55.5  15.0  0.00133   9.21 0.00178  
2 20937            54.8        59.5    49.3   5.57 0.00116   9.22 0.000212 
3 21109            44.6        58.5    48.8  -4.19 0.00125   9.22 0.000130 
4 21111            52.5        60.6    49.9   2.65 0.00107   9.22 0.0000445
5 21113            49.4        78.0    58.9  -9.48 0.00247   9.21 0.00131  
6 21114            48.1        49.4    44.0   4.15 0.00290   9.22 0.000295 
# ℹ 1 more variable: .std.resid <dbl>

augment()の出力

• .fitted: 予測値（\(\hat{y}\)）
• .resid: 残差（\(y - \hat{y}\)）
• .hat: てこ比（leverage）
• .cooksd: Cook’s distance（影響度）
• .std.resid: 標準化残差

augment()はpredict()と同様の機能を持つが，残差や診断指標も同時に得られる点で便利である．なお，augment(model, newdata = 新データ)で新しいデータに対する予測も可能．

# augment()を使った残差プロット
augment(model1) |>
  ggplot(aes(x = .fitted, y = .resid)) +
  geom_point(alpha = 0.3) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "予測値", y = "残差",
       title = "残差プロット（予測値 vs 残差）")

3.4.3 重回帰分析

複数の説明変数を同時に投入することで，他の変数を統制した上での関連を推定できる．

# 複数の説明変数
model2 <- lm(vocab_w2_g9 ~
               vocab_w2_g3 + sex + f_edu,
             data = d)
summary(model2)

Call:
lm(formula = vocab_w2_g9 ~ vocab_w2_g3 + sex + f_edu, data = d)

Residuals:
    Min      1Q  Median      3Q     Max 
-32.225  -5.657  -0.047   4.873  39.032 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)     18.60030    1.96726   9.455  < 2e-16 ***
vocab_w2_g3      0.50152    0.02991  16.768  < 2e-16 ***
sex男子         -0.66798    0.62034  -1.077 0.281867    
f_edu専門・短大 -1.03225    1.07219  -0.963 0.335939    
f_edu大学・院    2.42986    0.71893   3.380 0.000757 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.17 on 874 degrees of freedom
  (28967 個の観測値が欠損のため削除されました)
Multiple R-squared:  0.2762,    Adjusted R-squared:  0.2729 
F-statistic: 83.39 on 4 and 874 DF,  p-value: < 2.2e-16

tidy(model2)

# A tibble: 5 × 5
  term            estimate std.error statistic  p.value
  <chr>              <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)       18.6      1.97       9.45  2.89e-20
2 vocab_w2_g3        0.502    0.0299    16.8   6.29e-55
3 sex男子           -0.668    0.620     -1.08  2.82e- 1
4 f_edu専門・短大   -1.03     1.07      -0.963 3.36e- 1
5 f_edu大学・院      2.43     0.719      3.38  7.57e- 4

3.4.4 カテゴリカル変数の参照カテゴリ

カテゴリカル変数はダミー変数に変換されるが，どのカテゴリを参照（基準）とするかで係数の解釈が変わる．fct_relevel()で参照カテゴリを変更できる．

# 参照カテゴリを変更する場合
d <- d |>
  mutate(
    f_edu_f = fct_relevel(f_edu, "大学・院")
  )

model3 <- lm(vocab_w2_g9 ~ 
               vocab_w2_g3 + sex + f_edu_f, 
             data = d)
tidy(model3)

# A tibble: 5 × 5
  term              estimate std.error statistic  p.value
  <chr>                <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)         21.0      2.02       10.4  5.36e-24
2 vocab_w2_g3          0.502    0.0299     16.8  6.29e-55
3 sex男子             -0.668    0.620      -1.08 2.82e- 1
4 f_edu_f中学・高校   -2.43     0.719      -3.38 7.57e- 4
5 f_edu_f専門・短大   -3.46     0.986      -3.51 4.66e- 4

3.4.5 交互作用

説明変数の効果が他の変数によって異なる場合，交互作用項を投入する．*を使うと主効果と交互作用項の両方が含まれる．

# 性別×語彙力の交互作用
model4 <- lm(vocab_w2_g9 ~
               vocab_w2_g3 * sex,
             data = d)
tidy(model4)

# A tibble: 4 × 5
  term                estimate std.error statistic  p.value
  <chr>                  <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)          16.2       2.45        6.61 6.13e-11
2 vocab_w2_g3           0.555     0.0378     14.7  1.32e-44
3 sex男子               4.23      3.51        1.20 2.29e- 1
4 vocab_w2_g3:sex男子  -0.0691    0.0536     -1.29 1.98e- 1

3.4.6 回帰分析結果の図示

係数と信頼区間を可視化すると，効果の大きさや有意性を直感的に把握できる．

# 係数のプロット
tidy(model2, conf.int = TRUE) |>
  filter(term != "(Intercept)") |>
  ggplot(aes(x = estimate,
             y = term,
             xmin = conf.low,
             xmax = conf.high)) +
  geom_pointrange() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  labs(x = "係数（95%信頼区間）", y = "")

3.4.7 予測値の算出：predict()

回帰分析の結果から予測値を算出するpredict()関数は，様々な場面で活用される．

• 記述的な目的：モデルによる予測値と信頼区間を図示する
• 欠損値処理：回帰代入法で欠損値を予測値で補完する
• 因果推論：傾向スコアの算出，反実仮想の予測

3.4.7.1 基本的な使い方

# 分析用データの作成（欠損のないケース）
d_pred <- d |>
  drop_na(vocab_w2_g3, vocab_w2_g9, sex, f_edu)

# 回帰モデルの推定
model_pred <- lm(vocab_w2_g9 ~ vocab_w2_g3 + sex + f_edu,
                 data = d_pred)

# 予測値の算出
d_pred <- d_pred |>
  mutate(
    predicted = predict(model_pred),           # 予測値
    residual = vocab_w2_g9 - predicted         # 残差
  )

# 確認
d_pred |>
  select(vocab_w2_g9, predicted, residual) |>
  head()

# A tibble: 6 × 3
  vocab_w2_g9 predicted residual
        <dbl>     <dbl>    <dbl>
1        70.5      56.8    13.7 
2        54.8      50.9     3.96
3        44.6      46.2    -1.68
4        52.5      48.3     4.17
5        49.4      59.5   -10.1 
6        48.1      43.4     4.78

3.4.7.2 信頼区間付きの予測

predict()にinterval = "confidence"を指定すると，予測値の95%信頼区間が得られる．

# 信頼区間付きで予測
pred_ci <- predict(model_pred, interval = "confidence") |>
  as_tibble()

# データに結合
d_pred <- d_pred |>
  bind_cols(pred_ci |> rename(pred_fit = fit,
                               pred_lwr = lwr,
                               pred_upr = upr))

d_pred |>
  select(vocab_w2_g9, pred_fit, pred_lwr, pred_upr) |>
  head()

# A tibble: 6 × 4
  vocab_w2_g9 pred_fit pred_lwr pred_upr
        <dbl>    <dbl>    <dbl>    <dbl>
1        70.5     56.8     55.8     57.9
2        54.8     50.9     49.8     51.9
3        44.6     46.2     44.3     48.2
4        52.5     48.3     47.0     49.7
5        49.4     59.5     58.3     60.7
6        48.1     43.4     41.9     44.9

interval引数について

• "confidence": 平均予測値の信頼区間（回帰直線の不確実性）
• "prediction": 個々の観測値の予測区間（個人差も含む，より広い）

3.4.7.3 新しいデータに対する予測

newdata引数を使うと，元データにない値に対しても予測ができる．expand_grid()を使うと，指定した変数の全組み合わせを簡単に作成できる．

# 予測用の新しいデータを作成
# expand_grid()で全組み合わせを生成
new_data <- expand_grid(
  vocab_w2_g3 = seq(30, 70, by = 5),
  sex = c("男子", "女子"),
  f_edu = "大学・院"  # 固定する変数は1つの値だけ指定
)
new_data  # 9 × 2 = 18行のデータ

# A tibble: 18 × 3
   vocab_w2_g3 sex   f_edu   
         <dbl> <chr> <chr>   
 1          30 男子  大学・院
 2          30 女子  大学・院
 3          35 男子  大学・院
 4          35 女子  大学・院
 5          40 男子  大学・院
 6          40 女子  大学・院
 7          45 男子  大学・院
 8          45 女子  大学・院
 9          50 男子  大学・院
10          50 女子  大学・院
11          55 男子  大学・院
12          55 女子  大学・院
13          60 男子  大学・院
14          60 女子  大学・院
15          65 男子  大学・院
16          65 女子  大学・院
17          70 男子  大学・院
18          70 女子  大学・院

# 新しいデータに対する予測（信頼区間付き）
pred_new <- predict(model_pred,
                    newdata = new_data,
                    interval = "confidence") |>
  as_tibble()

# 結合
new_data <- new_data |>
  bind_cols(pred_new)

new_data

# A tibble: 18 × 6
   vocab_w2_g3 sex   f_edu      fit   lwr   upr
         <dbl> <chr> <chr>    <dbl> <dbl> <dbl>
 1          30 男子  大学・院  35.4  33.0  37.8
 2          30 女子  大学・院  36.1  33.8  38.4
 3          35 男子  大学・院  37.9  35.8  40.0
 4          35 女子  大学・院  38.6  36.5  40.6
 5          40 男子  大学・院  40.4  38.6  42.3
 6          40 女子  大学・院  41.1  39.3  42.9
 7          45 男子  大学・院  42.9  41.3  44.5
 8          45 女子  大学・院  43.6  42.1  45.1
 9          50 男子  大学・院  45.4  44.0  46.8
10          50 女子  大学・院  46.1  44.8  47.4
11          55 男子  大学・院  47.9  46.7  49.2
12          55 女子  大学・院  48.6  47.5  49.8
13          60 男子  大学・院  50.5  49.4  51.5
14          60 女子  大学・院  51.1  50.1  52.1
15          65 男子  大学・院  53.0  52.0  54.0
16          65 女子  大学・院  53.6  52.7  54.6
17          70 男子  大学・院  55.5  54.4  56.5
18          70 女子  大学・院  56.1  55.1  57.1

3.4.7.4 予測値の図示

# 予測値と信頼区間の図示
new_data |>
  ggplot(aes(x = vocab_w2_g3, y = fit,
             color = sex, fill = sex, linetype = sex)) +
  geom_ribbon(aes(ymin = lwr, ymax = upr),
              alpha = 0.2, color = NA) +
  geom_line() +
  scale_color_viridis_d(option = "C", end = .8) +
  scale_fill_viridis_d(option = "C", end = .8) +
  labs(x = "語彙力スコア（小3基準）",
       y = "読解力スコア（中3基準）の予測値",
       color = "性別", fill = "性別", linetype = "性別",
       title = "回帰モデルによる予測（父学歴：大学・院）") +
  theme(legend.position = "bottom")

predict()の応用

predict()は発展的な分析でも重要な役割を果たす．

• 回帰代入法（Regression Imputation）: 欠損値を回帰モデルの予測値で補完する，または因果推論で潜在アウトカム（反実仮想）を予測する
• 傾向スコア分析: ロジスティック回帰で処置確率（傾向スコア）を予測する

これらの手法では，推定したモデルを使って新しい条件下での値を予測することが本質的な役割となる．

3.5 ロジスティック回帰分析

2値（0/1）の従属変数を分析する場合は，ロジスティック回帰を使用する．

3.5.1 基本的なロジスティック回帰

glm()関数にfamily = binomialを指定する．係数は対数オッズ（log-odds）として推定される．

# 通塾（2015年時点）を従属変数とした分析
# 分析用データの作成
d_logit <- d |>
  filter(grade__2015 == 9) |>
  drop_na(sex, juku__2015, f_edu, income__2015) |>
  mutate(
    income_100 = income__2015 / 100  # 係数を見やすくするため100万円単位に
  )

# ロジスティック回帰
logit1 <- glm(juku__2015 ~ sex + f_edu + income_100,
              data = d_logit,
              family = binomial(link = "logit"))
summary(logit1)

Call:
glm(formula = juku__2015 ~ sex + f_edu + income_100, family = binomial(link = "logit"), 
    data = d_logit)

Coefficients:
                Estimate Std. Error z value Pr(>|z|)  
(Intercept)      0.08235    0.17524   0.470   0.6384  
sex男子          0.13082    0.12492   1.047   0.2950  
f_edu専門・短大  0.03919    0.19216   0.204   0.8384  
f_edu大学・院    0.03038    0.14692   0.207   0.8362  
income_100       0.04303    0.02012   2.139   0.0324 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1457.1  on 1095  degrees of freedom
Residual deviance: 1450.5  on 1091  degrees of freedom
AIC: 1460.5

Number of Fisher Scoring iterations: 4

3.5.2 結果の整理

tidy()のexponentiate = TRUEでオッズ比を取得できる．

# 係数（対数オッズ）
tidy(logit1, conf.int = TRUE, exponentiate = FALSE)

# A tibble: 5 × 7
  term            estimate std.error statistic p.value conf.low conf.high
  <chr>              <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 (Intercept)       0.0824    0.175      0.470  0.638  -0.262      0.425 
2 sex男子           0.131     0.125      1.05   0.295  -0.114      0.376 
3 f_edu専門・短大   0.0392    0.192      0.204  0.838  -0.336      0.418 
4 f_edu大学・院     0.0304    0.147      0.207  0.836  -0.258      0.318 
5 income_100        0.0430    0.0201     2.14   0.0324  0.00410    0.0831

# オッズ比として表示
tidy(logit1, conf.int = TRUE, exponentiate = TRUE)

# A tibble: 5 × 7
  term            estimate std.error statistic p.value conf.low conf.high
  <chr>              <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
1 (Intercept)         1.09    0.175      0.470  0.638     0.769      1.53
2 sex男子             1.14    0.125      1.05   0.295     0.892      1.46
3 f_edu専門・短大     1.04    0.192      0.204  0.838     0.715      1.52
4 f_edu大学・院       1.03    0.147      0.207  0.836     0.772      1.37
5 income_100          1.04    0.0201     2.14   0.0324    1.00       1.09

オッズ比の解釈

• オッズ比 > 1：その変数が増えると通塾確率が上がる
• オッズ比 < 1：その変数が増えると通塾確率が下がる
• オッズ比 = 1：関連なし

3.5.3 限界効果の推定

オッズ比は解釈が難しいため，限界効果（確率の変化）を計算すると分かりやすくなる．

# パッケージのインストール（初回のみ）
install.packages("marginaleffects")

library(marginaleffects)

# 平均限界効果（Average Marginal Effects: AME）
marginaleffects::avg_slopes(logit1)

       Term                Contrast Estimate Std. Error     z Pr(>|z|)   S
 f_edu      大学・院 - 中学・高校    0.00714    0.03455 0.207   0.8363 0.3
 f_edu      専門・短大 - 中学・高校  0.00920    0.04503 0.204   0.8381 0.3
 income_100 dY/dX                    0.01009    0.00468 2.155   0.0311 5.0
 sex        男子 - 女子              0.03069    0.02928 1.048   0.2945 1.8
     2.5 % 97.5 %
 -0.060582 0.0749
 -0.079051 0.0974
  0.000915 0.0193
 -0.026699 0.0881

Type: response

限界効果の解釈

• 連続変数：その変数が1単位増加したときの確率の変化（パーセントポイント）
• カテゴリ変数：参照カテゴリと比較したときの確率の差

3.5.4 限界効果の図示

# カテゴリ変数の限界効果をプロット
marginaleffects::avg_slopes(logit1,
                            variables = "f_edu") |>
  ggplot(aes(x = estimate, y = contrast,
             xmin = conf.low, xmax = conf.high)) +
  geom_pointrange() +
  geom_vline(xintercept = 0, linetype = "dashed") +
  labs(x = "通塾確率の差（パーセントポイント）",
       y = "父親の学歴（参照：中学・高校）",
       title = "父親学歴の限界効果")

# 収入の限界効果を収入水準別にプロット
plot_slopes(logit1, variables = "income_100", condition = "income_100") +
  labs(x = "世帯年収（100万円）",
       y = "限界効果（確率の変化）",
       title = "世帯年収の限界効果")

3.5.5 expand_grid()を使った予測

marginaleffectsパッケージのpredictions()関数とexpand_grid()を組み合わせると，特定の条件下での予測確率を柔軟に計算できる．

# 予測用データを作成
pred_data <- expand_grid(
  f_edu = c("中学・高校", "専門・短大", "大学・院"),
  income_100 = c(3, 5, 7, 10),  # 年収300万〜1000万
  sex = c("男子", "女子")
)

# 各条件での予測確率を計算
marginaleffects::predictions(logit1, newdata = pred_data) |>
  select(f_edu, income_100, sex, estimate, conf.low, conf.high) |>
  arrange(f_edu, income_100, sex)

 Estimate 2.5 % 97.5 %
    0.560 0.486  0.632
    0.592 0.519  0.662
    0.581 0.520  0.640
    0.613 0.552  0.670
    0.602 0.550  0.652
    0.633 0.581  0.682
    0.633 0.581  0.681
    0.662 0.610  0.711
    0.562 0.472  0.648
    0.594 0.507  0.676
    0.583 0.499  0.663
    0.615 0.533  0.691
    0.604 0.521  0.681
    0.635 0.555  0.709
    0.635 0.547  0.714
    0.664 0.578  0.741
    0.553 0.484  0.619
    0.585 0.518  0.648
    0.574 0.512  0.633
    0.605 0.545  0.663
    0.595 0.534  0.653
    0.626 0.566  0.682
    0.625 0.555  0.691
    0.656 0.587  0.718

expand_grid()の活用

expand_grid()で変数の全組み合わせを作成し，predictions()で予測確率を計算することで，「父学歴が大卒で年収700万の男子の通塾確率は何%か」といった具体的な問いに答えられる．これは因果推論における反実仮想の予測にも応用できる．

3.6 回帰分析結果の比較

複数のモデルを比較する際は，modelsummaryパッケージが便利である．

# パッケージのインストール（初回のみ）
install.packages("modelsummary")

library(modelsummary)

# 複数モデルの比較表
modelsummary(
  list("Model 1" = model1,
       "Model 2" = model2,
       "Model 4" = model4),
  stars = TRUE,
  gof_omit = "IC|Log|RMSE"
)

	Model 1	Model 2	Model 4
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
(Intercept)	18.321***	18.600***	16.197***
	(1.753)	(1.967)	(2.451)
vocab_w2_g3	0.520***	0.502***	0.555***
	(0.027)	(0.030)	(0.038)
sex男子		-0.668	4.231
		(0.620)	(3.512)
f_edu専門・短大		-1.032
		(1.072)
f_edu大学・院		2.430***
		(0.719)
vocab_w2_g3 × sex男子			-0.069
			(0.054)
Num.Obs.	1063	879	1062
R2	0.262	0.276	0.264
R2 Adj.	0.261	0.273	0.262

3.6.1 ロジスティック回帰の比較

ロジスティック回帰でも同様に比較表を作成できる．exponentiate = TRUEでオッズ比として表示する．

# モデルを追加
logit2 <- glm(juku__2015 ~ f_edu + income_100 + sex,
              data = d_logit,
              family = binomial(link = "logit"))

# 比較表（オッズ比で表示）
modelsummary(
  list("基本モデル" = logit1,
       "性別追加" = logit2),
  exponentiate = TRUE,
  stars = TRUE,
  gof_omit = "IC|Log"
)

	基本モデル	性別追加
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
(Intercept)	1.086	1.086
	(0.190)	(0.190)
sex男子	1.140	1.140
	(0.142)	(0.142)
f_edu専門・短大	1.040	1.040
	(0.200)	(0.200)
f_edu大学・院	1.031	1.031
	(0.151)	(0.151)
income_100	1.044*	1.044*
	(0.021)	(0.021)
Num.Obs.	1096	1096
F	1.602	1.602
RMSE	0.48	0.48

modelsummaryの利点

• 複数モデルを横に並べて比較できる
• HTML，LaTeX，Word形式で出力可能
• exponentiate = TRUEでオッズ比表示に切り替え

3.6.2 回帰分析における多重比較調整

回帰分析でカテゴリ変数を説明変数として使う場合，参照カテゴリとの差のみが出力される．すべてのカテゴリ間の比較を行い，多重比較調整をするにはemmeansパッケージが便利である．

# パッケージのインストール（初回のみ）
install.packages("emmeans")

library(emmeans)

# 父親学歴を説明変数とするモデル
model_edu <- lm(vocab_w2_g9 ~ f_edu, data = d)

# 推定周辺平均（Estimated Marginal Means）
emm <- emmeans(model_edu, ~ f_edu)
emm

      f_edu emmean    SE  df lower.CL upper.CL
 中学・高校   49.5 0.665 914     48.2     50.8
 専門・短大   48.5 1.010 914     46.5     50.4
 大学・院     53.5 0.445 914     52.7     54.4

Confidence level used: 0.95 

# すべてのペアワイズ比較（Holm法で多重比較調整）
pairs(emm, adjust = "holm")

           contrast      estimate   SE  df t.ratio p.value
 中学・高校 - 専門・短大     1.01 1.21 914   0.830  0.4068
 中学・高校 - 大学・院      -4.06 0.80 914  -5.074  <.0001
 専門・短大 - 大学・院      -5.07 1.11 914  -4.582  <.0001

P value adjustment: holm method for 3 tests 

emmeansの多重比較調整法

adjust引数で調整法を指定できる：

• "holm": Holm法（推奨．Bonferroniより検出力が高く，同じ第1種過誤率を維持）
• "tukey": Tukeyの方法（ANOVAの事後比較で伝統的）
• "fdr" または "BH": Benjamini-Hochberg法（多数の検定で偽発見率を制御）
• "bonferroni": Bonferroni法（保守的すぎることが多い）
• "none": 調整なし

推定周辺平均（EMM）とは

推定周辺平均は，他の共変量を平均値に固定したときの各カテゴリの予測平均である．複数の説明変数があるモデルでも，特定のカテゴリ変数の効果を分離して比較できる．

3.7 主成分分析と合成変数の活用

多数の変数がある場合，それらを少数の成分に要約する主成分分析（PCA: Principal Component Analysis）が有用である．相関の高い変数群を1つの成分にまとめることで，合成変数を作成できる．

ここでは学校適応に関する5項目から「学校適応得点」を作成し，グループ間の比較を行う例を示す．

3.7.1 主成分分析の実行

# 学校適応に関する項目を選択（Wave1子ども調査）
# ds0000000244_w1c: 授業が楽しい
# ds0000000245_w1c: 友だちとすごすのが楽しい
# ds0000000246_w1c: 尊敬できる先生がいる
# ds0000000249_w1c: 自分のクラスが好きだ
# ds0000000250_w1c: 自分の学校が好きだ
# 値: 1=とてもあてはまる, 2=まああてはまる, 3=あまりあてはまらない, 4=まったくあてはまらない

d_pca <- d |>
  select(
    school_1 = ds0000000244_w1c,
    school_2 = ds0000000245_w1c,
    school_3 = ds0000000246_w1c,
    school_4 = ds0000000249_w1c,
    school_5 = ds0000000250_w1c
  ) |>
  # 欠損値コード（7777, 8888, 9999）をNAに変換
  mutate(across(everything(), recode_missing)) |>
  drop_na()

# 主成分分析の実行
pca_result <- prcomp(d_pca, scale. = TRUE)

# 結果の要約（累積寄与率を確認）
summary(pca_result)

Importance of components:
                          PC1    PC2    PC3     PC4     PC5
Standard deviation     1.6992 0.8920 0.7363 0.70396 0.52847
Proportion of Variance 0.5775 0.1592 0.1084 0.09911 0.05586
Cumulative Proportion  0.5775 0.7366 0.8450 0.94414 1.00000

# 各主成分の負荷量（元の変数との関係）
pca_result$rotation

               PC1        PC2         PC3         PC4         PC5
school_1 0.4361184  0.4154288 -0.02879975 -0.79103168 -0.10324276
school_2 0.4040098 -0.5540224 -0.71889406 -0.05508858  0.09995907
school_3 0.3966815  0.6416295 -0.32029072  0.57194363 -0.03535093
school_4 0.4888489 -0.3050556  0.39995292  0.18459877 -0.68842489
school_5 0.5002472 -0.1254223  0.46884264  0.10018937  0.71004972

scale.引数について

scale. = TRUEを指定すると，各変数を標準化（平均0，標準偏差1）してから分析する．変数間でスケールが異なる場合は標準化が推奨される．

3.7.2 主成分得点をデータに追加

# 主成分得点をデータに追加
# 欠損値がある行はNAになる
d_school <- d |>
  select(
    school_1 = ds0000000244_w1c,
    school_2 = ds0000000245_w1c,
    school_3 = ds0000000246_w1c,
    school_4 = ds0000000249_w1c,
    school_5 = ds0000000250_w1c
  ) |>
  mutate(across(everything(), recode_missing))

# 主成分得点を元のデータに追加（欠損がある行は自動的にNAになる）
d <- d |>
  mutate(school_adapt = predict(pca_result, newdata = d_school)[, 1])

# 学校適応得点の分布
d |>
  drop_na(school_adapt) |>
  ggplot(aes(x = school_adapt)) +
  geom_histogram(binwidth = 0.5, fill = "steelblue", color = "white") +
  labs(x = "学校適応得点（第1主成分）", y = "度数",
       title = "学校適応得点の分布")

3.7.3 合成変数を用いたグループ間比較

作成した学校適応得点を使って，父親学歴別の平均値と95%信頼区間を算出・図示する．

# 父親学歴別の学校適応得点（平均と95%信頼区間）
d |>
  drop_na(school_adapt, f_edu) |>
  summarise(
    mean = mean(school_adapt),
    sd = sd(school_adapt),
    n = n(),
    se = sd / sqrt(n),
    ll = mean - 1.96 * se,
    ul = mean + 1.96 * se,
    .by = f_edu
  )

# A tibble: 3 × 7
  f_edu         mean    sd     n     se      ll      ul
  <fct>        <dbl> <dbl> <int>  <dbl>   <dbl>   <dbl>
1 大学・院   -0.105   1.66  5590 0.0222 -0.148  -0.0614
2 中学・高校  0.0913  1.72  3323 0.0298  0.0328  0.150 
3 専門・短大 -0.0217  1.73  1572 0.0436 -0.107   0.0637

# 図示（信頼区間付き）
d |>
  drop_na(school_adapt, f_edu) |>
  summarise(
    mean = mean(school_adapt),
    sd = sd(school_adapt),
    n = n(),
    se = sd / sqrt(n),
    ll = mean - 1.96 * se,
    ul = mean + 1.96 * se,
    .by = f_edu
  ) |>
  ggplot(aes(x = f_edu, y = mean, ymin = ll, ymax = ul)) +
  geom_col(fill = "steelblue", alpha = 0.7) +
  geom_errorbar(width = 0.2) +
  labs(x = "父親の学歴", y = "学校適応得点（平均）",
       title = "父親学歴別の学校適応得点（95%信頼区間付き）")

3.7.4 合成変数を用いた回帰分析

学校適応得点を説明変数として回帰分析に使用することもできる．

# 学校適応得点が通塾に与える影響（ロジスティック回帰）
model_adapt <- glm(juku__2015 ~ school_adapt + sex + f_edu,
                   data = d,
                   family = binomial)

tidy(model_adapt, conf.int = TRUE) |>
  mutate(
    odds_ratio = exp(estimate),
    or_ll = exp(conf.low),
    or_ul = exp(conf.high)
  ) |>
  select(term, odds_ratio, or_ll, or_ul, p.value)

# A tibble: 5 × 5
  term            odds_ratio or_ll or_ul   p.value
  <chr>                <dbl> <dbl> <dbl>     <dbl>
1 (Intercept)          0.395 0.363 0.430 4.87e-101
2 school_adapt         1.00  0.979 1.03  8.24e-  1
3 sex男子              1.09  1.00  1.18  4.50e-  2
4 f_edu専門・短大      1.30  1.14  1.48  5.65e-  5
5 f_edu大学・院        1.77  1.61  1.94  1.00e- 33

主成分分析の活用場面

• 次元削減：多数の変数を少数の主成分に要約し，回帰分析の説明変数として使用
• 多重共線性の回避：相関の高い変数群を1つの成分にまとめる
• 尺度得点の作成：複数の質問項目から合成得点を作成

関連手法

• 探索的因子分析（EFA）：観測変数の背後にある潜在因子を探索する手法．主成分分析と異なり、観測変数の共通因子と独自因子を区別する．psychパッケージのfa()関数が推奨（oblimin等の斜交回転に対応）．base Rのfactanal()はvarimax/promaxのみ
• 対応分析（CA）：カテゴリ変数のクロス表を可視化する手法．パッケージ：FactoMineR, ca
• 多重対応分析（MCA）：複数のカテゴリ変数を同時に分析する手法．パッケージ：FactoMineR, homals

主成分分析と因子分析の違い

• 主成分分析（PCA）：データの分散を最大限説明する成分を抽出。次元削減・合成変数作成が目的。
• 因子分析（FA）：観測変数の共分散構造を説明する潜在因子を推定。構成概念の探索が目的。

実用上、項目数が多く因子構造の探索が目的なら因子分析、単純な次元削減・合成変数作成なら主成分分析を使うことが多い。

3.8 ウェイト付き分析

調査データを分析する際，サンプリングウェイト（重み付け）を用いて母集団への一般化を行う場合がある。

3.8.1 weights引数を使う方法

多くの関数にはweights引数があり，単純なウェイト付き分析が可能である。

# 仮のウェイト変数を作成（実際のデータではウェイト変数を使用）
d <- d |>
  mutate(weight = runif(n(), 0.5, 1.5))

# ウェイト付き回帰分析
model_weighted <- lm(vocab_w2_g9 ~ f_edu + income__2015,
                     data = d,
                     weights = weight)
summary(model_weighted)

# ウェイト付きロジスティック回帰
model_logit_weighted <- glm(juku__2015 ~ f_edu + income__2015,
                            data = d,
                            family = binomial,
                            weights = weight)
summary(model_logit_weighted)

weights引数の限界

lm()やglm()のweights引数は、標準誤差の計算において複雑なサンプリングデザイン（層化、クラスタリングなど）を考慮しない。単純無作為抽出の場合や、重み付けした点推定のみが必要な場合に使用する。

3.8.2 surveyパッケージを使う方法

複雑なサンプリングデザイン（層化抽出、クラスタ抽出など）を考慮した分析にはsurveyパッケージを使用する。

library(survey)

# サンプリングデザインの指定
# ids: クラスタ変数（なければ~1）
# strata: 層化変数（なければNULL）
# weights: ウェイト変数
# data: データフレーム
design <- svydesign(
  ids = ~1,           # クラスタなし（単純無作為抽出）
  strata = NULL,      # 層化なし
  weights = ~weight,  # ウェイト変数
  data = d
)

# ウェイト付き平均
svymean(~vocab_w2_g9, design, na.rm = TRUE)

# ウェイト付き集計（カテゴリ別）
svyby(~vocab_w2_g9, ~f_edu, design, svymean, na.rm = TRUE)

# ウェイト付きクロス表
svytable(~f_edu + juku__2015, design)

# ウェイト付き回帰分析
model_survey <- svyglm(vocab_w2_g9 ~ f_edu + income__2015, design)
summary(model_survey)

# ウェイト付きロジスティック回帰
model_survey_logit <- svyglm(juku__2015 ~ f_edu + income__2015,
                              design,
                              family = quasibinomial)
summary(model_survey_logit)

svyglm()でのロジスティック回帰

svyglm()でロジスティック回帰を行う場合は，family = quasibinomialを指定する。binomialを使うと収束の問題が生じる場合がある。

3.8.3 層化・クラスタ抽出の場合

# 層化クラスタ抽出の例
# 例：都道府県で層化し、市区町村でクラスタ抽出
design_complex <- svydesign(
  ids = ~cluster_id,      # クラスタ変数（市区町村など）
  strata = ~stratum_id,   # 層化変数（都道府県など）
  weights = ~weight,
  data = d,
  nest = TRUE             # クラスタが層内でネストされている場合
)

# 分析は同様
svymean(~vocab_w2_g9, design_complex, na.rm = TRUE)

srvyrパッケージ

srvyrパッケージを使うと、surveyパッケージの機能をtidyverseの文法で使える。

library(srvyr)

# srvyrでデザイン指定
design_srvyr <- d |>
  as_survey_design(weights = weight)

# dplyr風の集計
design_srvyr |>
  group_by(f_edu) |>
  summarise(
    mean_vocab = survey_mean(vocab_w2_g9, na.rm = TRUE),
    n = unweighted(n())
  )

ウェイトを使うべきかどうか

ウェイト付き分析が常に必要なわけではない。

• 記述統計（母集団の特性を推定）：ウェイトを使用
• 回帰分析（変数間の関連の推定）：ウェイトが必要かは議論がある

回帰分析でウェイトを使うかどうかは、分析の目的やサンプリングデザインに依存する。Solon, Haider, & Wooldridge (2015)などの議論を参照のこと。

3.9 パネルデータの作成

パネルデータ（縦断データ）では，同一個人を複数時点で観察する．

3.9.1 ワイド形式からロング形式への変換

パネルデータには2つの形式がある．

• ワイド形式：1行に1人のデータ，各時点の変数が横に並ぶ
• ロング形式：1行に1人1時点のデータ，時点ごとに行が増える

3.9.1.1 ワイド形式（Wide format）

1行に1人のすべての時点のデータが含まれる．

| id | income__2015 | income__2016 | income__2017 |
|----|--------------|--------------|--------------|
| 1  | 500          | 550          | 600          |
| 2  | 400          | 420          | 450          |
| 3  | 700          | 710          | 750          |

3.9.1.2 ロング形式（Long format）

1行に1人1時点のデータが含まれ，時点ごとに行が増える．

| id | year | income |
|----|------|--------|
| 1  | 2015 | 500    |
| 1  | 2016 | 550    |
| 1  | 2017 | 600    |
| 2  | 2015 | 400    |
| 2  | 2016 | 420    |
| 2  | 2017 | 450    |
| 3  | 2015 | 700    |
| 3  | 2016 | 710    |
| 3  | 2017 | 750    |

多くの統計手法（固定効果モデルなど）ではロング形式が必要である．pivot_longer()関数で変換できる．

変数名の命名規則

ワイド形式からロング形式に変換する際，変数名を変数名__年のように統一しておくと，names_sep = "__"で簡単に分割できる．

d_wide <- d |>
  select(PanelID,
         sex,
         f_edu,
         m_edu,
         contains("grade__"),
         contains("income__"),
         contains("juku__"))

d_long <- d_wide |>
    pivot_longer(
    cols = -c(PanelID, sex, f_edu, m_edu),
    names_to = c(".value", "year"),
    names_sep = "__"
  ) |>
  mutate(year = as.numeric(year))

# 確認
d_long |> head(10)

# A tibble: 10 × 8
   PanelID sex   f_edu      m_edu       year grade income  juku
     <dbl> <chr> <fct>      <fct>      <dbl> <dbl>  <dbl> <dbl>
 1 1600001 男子  専門・短大 専門・短大  2015     0     NA    NA
 2 1600001 男子  専門・短大 専門・短大  2016     1    700     0
 3 1600001 男子  専門・短大 専門・短大  2017     2    700     0
 4 1600001 男子  専門・短大 専門・短大  2018     3    700     0
 5 1600001 男子  専門・短大 専門・短大  2019     4    700     0
 6 1600001 男子  専門・短大 専門・短大  2020     5    550     0
 7 1600001 男子  専門・短大 専門・短大  2021     6    700     1
 8 1600002 男子  <NA>       専門・短大  2015     0     NA    NA
 9 1600002 男子  <NA>       専門・短大  2016     1    100     0
10 1600002 男子  <NA>       専門・短大  2017     2    100     0

3.9.2 ロング形式でのデータ操作

ロング形式では，.by引数を使って個人別や年別の集計が容易にできる．

# 個人別の平均収入
d_long |>
  summarise(mean_income = mean(income, na.rm = TRUE),
            .by = PanelID) |>
  drop_na() |>
  head(10)

# A tibble: 10 × 2
   PanelID mean_income
     <dbl>       <dbl>
 1 1600001        675 
 2 1600002        100 
 3 1600003        690 
 4 1600004       1040 
 5 1600005        467.
 6 1600006        500 
 7 1600007        558.
 8 1600008       1250 
 9 1600009        608.
10 1600010        150 

# 年別の平均収入
d_long |>
  summarise(mean_income = mean(income, na.rm = TRUE),
            n = sum(!is.na(income)),
            .by = year) |>
  arrange(year)

# A tibble: 7 × 3
   year mean_income     n
  <dbl>       <dbl> <int>
1  2015        689. 14959
2  2016        704. 14353
3  2017        713. 13798
4  2018        729. 13147
5  2019        742. 13615
6  2020        748. 13790
7  2021        761. 13734

# 年別の通塾率
d_long |>
  summarise(juku_rate = mean(juku, na.rm = TRUE),
            n = sum(!is.na(juku)),
            .by = year) |>
  arrange(year)

# A tibble: 7 × 3
   year juku_rate     n
  <dbl>     <dbl> <int>
1  2015     0.312 16708
2  2016     0.300 15769
3  2017     0.310 15346
4  2018     0.300 14255
5  2019     0.294 15359
6  2020     0.294 15597
7  2021     0.307 15688

3.9.2.1 ラグ変数・リード変数の作成

パネルデータ分析では，前時点の値（ラグ変数）や次時点の値（リード変数）を使うことが多い．dplyrのlag()とlead()関数で作成できる．

# ラグ変数とリード変数の作成
d_long <- d_long |>
  arrange(PanelID, year) |>  # 必ずソートしてから

  mutate(
    income_lag = lag(income),    # 前年の収入
    income_lead = lead(income),  # 翌年の収入
    income_diff = income - income_lag,  # 収入の変化
    .by = PanelID  # 個人ごとに計算
  )

# 確認
d_long |>
  filter(PanelID == first(PanelID)) |>
  select(PanelID, year, income, income_lag, income_lead, income_diff)

# A tibble: 7 × 6
  PanelID  year income income_lag income_lead income_diff
    <dbl> <dbl>  <dbl>      <dbl>       <dbl>       <dbl>
1 1400001  2015    900         NA          NA          NA
2 1400001  2016     NA        900          NA          NA
3 1400001  2017     NA         NA          NA          NA
4 1400001  2018     NA         NA          NA          NA
5 1400001  2019     NA         NA          NA          NA
6 1400001  2020     NA         NA          NA          NA
7 1400001  2021     NA         NA          NA          NA

lag()使用時の注意

• 必ず.by（またはgroup_by()）で個人を指定：指定しないと他の個人のデータが混入する
• 事前にソート：arrange()で個人・時点の順に並べてから実行
• 欠損値の発生：最初の時点のラグ変数，最後の時点のリード変数はNAになる

3.9.3 パネルデータの可視化

パネルデータでは，時間変化をグラフで示すことが重要である．まず全体の推移を確認し，次にグループ別の違いを検討する．

# 年別の通塾率の推移（信頼区間付き）
d_long |>
  drop_na(juku) |>
  summarise(
    juku_rate = mean(juku, na.rm = TRUE),
    n = sum(!is.na(juku)),
    se = sqrt(juku_rate * (1 - juku_rate) / n),  # 二項分布の標準誤差
    ll = juku_rate - 1.96 * se,
    ul = juku_rate + 1.96 * se,
    .by = year) |>
  ggplot(aes(x = year, y = juku_rate, ymin = ll, ymax = ul)) +
  geom_line() +
  geom_point() +
  geom_ribbon(alpha = 0.2) +
  scale_x_continuous(breaks = 2015:2021) +
  scale_y_continuous(labels = scales::percent) +
  labs(x = "年", y = "通塾率",
       title = "通塾率の推移（95%信頼区間付き）")

# 父親学歴別の通塾率の推移（信頼区間付き）
d_long |>
  drop_na(juku, f_edu) |>
  summarise(
    juku_rate = mean(juku, na.rm = TRUE),
    n = sum(!is.na(juku)),
    se = sqrt(juku_rate * (1 - juku_rate) / n),
    ll = juku_rate - 1.96 * se,
    ul = juku_rate + 1.96 * se,
    .by = c(year, f_edu)) |>
  ggplot(aes(x = year, y = juku_rate, color = f_edu, fill = f_edu,
             linetype = f_edu, shape = f_edu,
             ymin = ll, ymax = ul)) +
  geom_ribbon(alpha = 0.15, color = NA) +
  geom_line() +
  geom_point() +
  scale_x_continuous(breaks = 2015:2021) +
  scale_y_continuous(labels = scales::percent) +
  scale_color_viridis_d(option = "C", end = .8) +
  scale_fill_viridis_d(option = "C", end = .8) +
  labs(x = "年", y = "通塾率",
       color = "父親の学歴", fill = "父親の学歴",
       linetype = "父親の学歴", shape = "父親の学歴",
       title = "父親学歴別の通塾率の推移（95%信頼区間付き）")

# 個人の収入変化（サンプル100人）
set.seed(123)
sample_ids <- d_long |>
  drop_na(income) |>
  distinct(PanelID) |>
  slice_sample(n = 100) |>
  pull(PanelID)

d_long |>
  filter(PanelID %in% sample_ids) |>
  ggplot(aes(x = year,
             y = income,
             group = PanelID)) +
  geom_line(alpha = 0.3) +
  geom_point(alpha = 0.5, size = 1) +
  stat_summary(aes(group = 1), fun = mean, geom = "line",
               color = "red", linewidth = 1.5) +
  stat_summary(aes(group = 1), fun = mean, geom = "point",
               color = "red", size = 3) +
  scale_x_continuous(breaks = 2015:2021) +
  labs(x = "年", y = "世帯年収（万円）",
       title = "世帯年収の推移（サンプル100人）",
       caption = "赤線は全体の平均値")

3.9.3.1 状態変化の可視化：Sankey図

パネルデータでは，個人の状態変化（例：通塾の有無）を追跡できる．ggalluvialパッケージを使うと，状態間の遷移をSankey図（alluvial plot）として可視化できる．

ここでは，中学1年から中学3年までの3年間で通塾状態がどう変化するかを可視化する．W1で小6だった児童を中1〜中3まで追跡する．

library(ggalluvial)

# W1で小6（2016年度に中1）の児童を中1〜中3まで追跡
# grade = 6 in W1 → grade 7 in W2 (2016), grade 8 in W3 (2017), grade 9 in W4 (2018)
d_sankey <- d |>
  filter(`w1学年` == 6) |>  # W1で小6

  select(PanelID, juku__2016, juku__2017, juku__2018) |>
  drop_na() |>
  mutate(
    中1 = if_else(juku__2016 == 1, "通塾", "非通塾"),
    中2 = if_else(juku__2017 == 1, "通塾", "非通塾"),
    中3 = if_else(juku__2018 == 1, "通塾", "非通塾")
  ) |>
  count(中1, 中2, 中3)

d_sankey

# A tibble: 8 × 4
  中1    中2    中3        n
  <chr>  <chr>  <chr>  <int>
1 通塾   通塾   通塾     243
2 通塾   通塾   非通塾    20
3 通塾   非通塾 通塾      16
4 通塾   非通塾 非通塾    15
5 非通塾 通塾   通塾     114
6 非通塾 通塾   非通塾    17
7 非通塾 非通塾 通塾     111
8 非通塾 非通塾 非通塾   291

# Sankey図
d_sankey |>
  ggplot(aes(axis1 = 中1, axis2 = 中2, axis3 = 中3, y = n)) +
  geom_alluvium(aes(fill = 中1), width = 1/12) +
  geom_stratum(width = 1/12, fill = "grey80", color = "grey50") +
  geom_text(stat = "stratum", aes(label = after_stat(stratum))) +
  scale_x_discrete(limits = c("中1", "中2", "中3"), expand = c(0.1, 0.1)) +
  scale_fill_manual(values = c("通塾" = "#0d0887", "非通塾" = "#f0f921")) +
  labs(y = "人数", fill = "中1時点の状態",
       title = "中学校3年間の通塾状態の変化") +
  theme(legend.position = "bottom")

Sankey図の活用場面

• 職業移動や階層移動の可視化
• 学歴・進路の変化
• 就業状態（就業・失業・非労働力）の遷移
• 健康状態や生活習慣の変化

3.9.4 学年を横軸にした可視化

通塾のように学年によって大きく変化する現象は，調査年ではなく学年を横軸にした方が理解しやすい．同じ学年でも調査年が異なるコーホートが存在するため，コーホート間の比較も可能になる．

# 学年別の通塾率
d_long |>
  drop_na(juku, grade) |>
  filter(grade %in% 1:12) |>
  summarise(juku_rate = mean(juku, na.rm = TRUE),
            n = sum(!is.na(juku)),
            .by = grade) |>
  arrange(grade)

# A tibble: 12 × 3
   grade juku_rate     n
   <dbl>     <dbl> <int>
 1     1     0.145 12451
 2     2     0.171 11534
 3     3     0.199 10895
 4     4     0.254  9842
 5     5     0.314  9349
 6     6     0.381  8839
 7     7     0.394  8157
 8     8     0.476  7931
 9     9     0.607  7940
10    10     0.222  7317
11    11     0.265  7185
12    12     0.385  7268

# 学年別の通塾率（図示・信頼区間付き）
d_long |>
  drop_na(juku, grade) |>
  filter(grade %in% 1:12) |>
  summarise(
    juku_rate = mean(juku, na.rm = TRUE),
    n = sum(!is.na(juku)),
    se = sqrt(juku_rate * (1 - juku_rate) / n),
    ll = juku_rate - 1.96 * se,
    ul = juku_rate + 1.96 * se,
    .by = grade) |>
  ggplot(aes(x = grade, y = juku_rate, ymin = ll, ymax = ul)) +
  geom_ribbon(alpha = 0.2) +
  geom_line() +
  geom_point() +
  scale_x_continuous(breaks = 1:12,
                     labels = c(paste0("小", 1:6), paste0("中", 1:3), paste0("高", 1:3))) +
  scale_y_continuous(labels = scales::percent) +
  labs(x = "学年", y = "通塾率",
       title = "学年別の通塾率（95%信頼区間付き）")

# 出生コーホート別の通塾率（同一コーホートの学年変化を追跡）
# コーホート = 調査年 - 学年（小1入学年度で定義）
d_long |>
  drop_na(juku, grade) |>
  filter(grade %in% 1:12) |>
  mutate(cohort = year - grade) |>
  summarise(
    juku_rate = mean(juku, na.rm = TRUE),
    n = sum(!is.na(juku)),
    .by = c(grade, cohort)) |>
  filter(n >= 30) |>  # サンプルサイズが小さいコーホート×学年を除外
  ggplot(aes(x = grade, y = juku_rate,
             color = cohort, group = cohort)) +
  geom_line() +
  geom_point() +
  scale_x_continuous(breaks = 1:12,
                     labels = c(paste0("小", 1:6), paste0("中", 1:3), paste0("高", 1:3))) +
  scale_y_continuous(labels = scales::percent) +
  scale_color_viridis_c(option = "C", end = 0.8) +
  labs(x = "学年", y = "通塾率",
       color = "出生コーホート\n（小1入学年度）",
       title = "出生コーホート別の通塾率の変化")

コーホート別分析の意義

調査年別の比較では異なる世代を比較しているが，出生コーホート別の分析では同じ個人集団の学年による変化を追跡できる．これはパネルデータの強みを活かした分析である．各線は同一コーホートの成長に伴う通塾率の変化を示している．

# 学年別・父親学歴別の通塾率（信頼区間付き）
d_long |>
  drop_na(juku, grade, f_edu) |>
  filter(grade %in% 1:12) |>
  summarise(
    juku_rate = mean(juku, na.rm = TRUE),
    n = sum(!is.na(juku)),
    se = sqrt(juku_rate * (1 - juku_rate) / n),
    ll = juku_rate - 1.96 * se,
    ul = juku_rate + 1.96 * se,
    .by = c(grade, f_edu)) |>
  ggplot(aes(x = grade, y = juku_rate,
             color = f_edu, fill = f_edu, linetype = f_edu,
             shape = f_edu, group = f_edu,
             ymin = ll, ymax = ul)) +
  geom_ribbon(alpha = 0.15, color = NA) +
  geom_line() +
  geom_point() +
  scale_x_continuous(breaks = 1:12,
                     labels = c(paste0("小", 1:6), paste0("中", 1:3), paste0("高", 1:3))) +
  scale_y_continuous(labels = scales::percent) +
  scale_color_viridis_d(option = "C", end = 0.8) +
  scale_fill_viridis_d(option = "C", end = 0.8) +
  labs(x = "学年", y = "通塾率",
       color = "父親の学歴", fill = "父親の学歴",
       linetype = "父親の学歴", shape = "父親の学歴",
       title = "学年別・父親学歴別の通塾率（95%信頼区間付き）")

横軸の選択：調査年 vs 学年

• 調査年を横軸：時代効果（コロナ禍など）を見たいとき，経済指標など年で変化する変数
• 学年を横軸：発達や教育段階による変化を見たいとき，通塾・進学希望など学年で変化する変数

調査年と学年の両方を使う場合，一方を横軸，他方を色や形（color, shape）で表現する．同じ学年でも調査年が異なるデータがあれば，コーホート効果の検討が可能になる．

3.9.5 map()を使った学年別分析

purrrパッケージのmap()関数を使うと，グループごとに同じ分析を繰り返し適用できる．ここでは学年別に重回帰分析を行い，係数を比較する．

# 学年別に重回帰分析を実行
grade_models <- d_long |>
  filter(grade %in% 1:12) |>
  drop_na(juku, income, f_edu, sex) |>
  nest(.by = grade) |>
  mutate(
    model = map(data, ~ glm(juku ~ income + f_edu + sex,
                            data = .x, family = binomial)),
    result = map(model, ~ tidy(.x, conf.int = TRUE))  # 信頼区間も取得
  ) |>
  unnest(result) |>
  select(grade, term, estimate, std.error, conf.low, conf.high, p.value)

grade_models

# A tibble: 60 × 7
   grade term             estimate std.error  conf.low conf.high   p.value
   <dbl> <chr>               <dbl>     <dbl>     <dbl>     <dbl>     <dbl>
 1    12 (Intercept)     -1.81     0.0860    -1.98      -1.64    3.17e- 98
 2    12 income           0.000988 0.0000821  0.000828   0.00115 2.36e- 33
 3    12 f_edu専門・短大  0.245    0.0959     0.0563     0.433   1.07e-  2
 4    12 f_edu大学・院    0.833    0.0701     0.696      0.971   1.35e- 32
 5    12 sex男子          0.125    0.0582     0.0113     0.240   3.13e-  2
 6    10 (Intercept)     -2.29     0.0948    -2.47      -2.10    2.47e-128
 7    10 income           0.000917 0.0000869  0.000747   0.00109 5.09e- 26
 8    10 f_edu専門・短大  0.166    0.111     -0.0535     0.383   1.36e-  1
 9    10 f_edu大学・院    0.515    0.0822     0.355      0.678   3.68e- 10
10    10 sex男子          0.00850  0.0658    -0.121      0.138   8.97e-  1
# ℹ 50 more rows

# 収入の係数を学年別にプロット（95%信頼区間付き）
grade_models |>
  filter(term == "income") |>
  ggplot(aes(x = grade, y = estimate,
             ymin = conf.low, ymax = conf.high)) +
  geom_ribbon(alpha = 0.2) +
  geom_line() +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  scale_x_continuous(breaks = 1:12,
                     labels = c(paste0("小", 1:6), paste0("中", 1:3), paste0("高", 1:3))) +
  labs(x = "学年", y = "収入の係数（対数オッズ）",
       title = "学年別：収入と通塾の関連（95%信頼区間付き）")

対数オッズ係数のモデル間比較の注意点

ロジスティック回帰の対数オッズ係数は，厳密にはモデル間で直接比較できない．これは，ロジスティック回帰では残差分散が識別できず，各モデルで異なるスケールになるためである（Mood, 2010）．

グループ間（ここでは学年間）で効果の大きさを比較したい場合は，平均限界効果（Average Marginal Effects: AME）を用いることが推奨される．AMEは確率スケールでの効果を表すため，モデル間で直接比較可能である．

# 平均限界効果（AME）を計算
# marginaleffectsパッケージを使用
library(marginaleffects)

# 学年別にAMEを計算
grade_ame <- d_long |>
  filter(grade %in% 1:12) |>
  drop_na(juku, income, f_edu, sex) |>
  nest(.by = grade) |>
  mutate(
    model = map(data, ~ glm(juku ~ income + f_edu + sex,
                            data = .x, family = binomial)),
    # 収入のAMEを計算
    ame = map(model, ~ avg_slopes(.x, variables = "income"))
  ) |>
  unnest(ame) |>
  select(grade, term, estimate, std.error, conf.low, conf.high)

grade_ame

# A tibble: 12 × 6
   grade term    estimate  std.error  conf.low conf.high
   <dbl> <chr>      <dbl>      <dbl>     <dbl>     <dbl>
 1    12 income 0.000216  0.0000171  0.000183   0.000250
 2    10 income 0.000155  0.0000142  0.000127   0.000183
 3     7 income 0.000150  0.0000174  0.000116   0.000184
 4     6 income 0.000284  0.0000165  0.000252   0.000316
 5     4 income 0.000231  0.0000129  0.000206   0.000256
 6    11 income 0.000160  0.0000155  0.000130   0.000190
 7     5 income 0.000268  0.0000144  0.000240   0.000296
 8     9 income 0.0000916 0.0000187  0.0000549  0.000128
 9     8 income 0.000138  0.0000186  0.000102   0.000175
10     3 income 0.000169  0.0000110  0.000147   0.000190
11     2 income 0.000145  0.0000100  0.000125   0.000164
12     1 income 0.000128  0.00000878 0.000111   0.000145

# AMEを学年別にプロット
grade_ame |>
  ggplot(aes(x = grade, y = estimate,
             ymin = conf.low, ymax = conf.high)) +
  geom_ribbon(alpha = 0.2) +
  geom_line() +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  scale_x_continuous(breaks = 1:12,
                     labels = c(paste0("小", 1:6), paste0("中", 1:3), paste0("高", 1:3))) +
  labs(x = "学年", y = "収入の限界効果（確率の変化）",
       title = "学年別：収入と通塾の関連（平均限界効果）",
       subtitle = "収入が1単位増加したときの通塾確率の変化")

map()の利点

map()を使うことで：

• 同じ分析を複数のグループに一括適用できる
• 結果を整理されたデータフレームとして取得できる
• グループ間の係数の比較が容易になる

nest()でグループごとにデータを分割し，map()で各グループにモデルを適用，unnest()で結果を展開するのが基本的な流れである．

3.10 パネルデータの分析

3.10.1 固定効果モデルの考え方

個人の観察されない異質性（ability, motivation など）を統制するために，固定効果モデルを使用する．

個人固定効果モデル（one-way）

\[ Y_{it} = \alpha_i + \beta X_{it} + \epsilon_{it} \]

• \(\alpha_i\)：個人固定効果（時間によって変化しない個人の特性）
• \(\beta\)：関心のあるパラメータ

個人＋時間固定効果モデル（two-way）

\[ Y_{it} = \alpha_i + \lambda_t + \beta X_{it} + \epsilon_{it} \]

• \(\alpha_i\)：個人固定効果（時間によって変化しない個人の特性）
• \(\lambda_t\)：時間固定効果（個人によらない時点ごとの影響，例：学習指導要領の改訂，コロナ禍による休校）
• \(\beta\)：関心のあるパラメータ

two-way固定効果モデルでは，個人に共通する時間トレンド（例：教育政策の変更やコロナ禍の影響）も統制できる．

3.10.2 plmパッケージによる固定効果モデル

plmパッケージは伝統的なパネルデータ分析用のパッケージである．pdata.frame()でパネル構造を指定し，model = "within"で固定効果モデルを推定する．

# install.packages("plm")
library(plm)

# パネルデータとして設定
pdata <- pdata.frame(d_long,
                     index = c("PanelID", "year"))

# 固定効果モデル（通塾と収入の関連）
fe_model <- plm(income ~ juku,
                data = pdata,
                model = "within")
summary(fe_model)

3.10.3 fixestパッケージによる固定効果モデル

fixestパッケージは高速で，多くの固定効果を扱える．

# install.packages("fixest")
library(fixest)

# 個人固定効果モデル
fe_model <- feols(income ~ juku | PanelID, data = d_long)
summary(fe_model)

# 個人固定効果と年固定効果
fe_model2 <- feols(income ~ juku | PanelID + year, data = d_long)
summary(fe_model2)

# 結果の比較
etable(fe_model, fe_model2)

fixestの利点

• 計算が非常に高速（大規模データに最適）
• 複数の固定効果を簡単に指定できる（|の後に+で追加）
• etable()で複数モデルの結果を比較できる
• クラスタロバスト標準誤差が標準で計算される

3.11 データの保存・エクスポート

加工したデータや分析結果を保存する方法を紹介する．

3.11.1 データフレームの保存

用途に応じてファイル形式を選択する．

# CSVファイルとして保存
write_csv(d_long, here("data", "processed", "d_long.csv"))

# Stataファイル（.dta）として保存
haven::write_dta(d_long, here("data", "processed", "d_long.dta"))

# RDS形式で保存（R専用，高速・圧縮）
saveRDS(d_long, here("data", "processed", "d_long.rds"))

# RDS形式で読み込む場合
# d_long <- readRDS(here("data", "processed", "d_long.rds"))

ファイル形式の選び方

形式	利点	欠点
CSV	汎用性が高い，他ソフトで開ける	ラベル情報が失われる
DTA	Stataで読める，ラベル保持	Stata専用
RDS	高速，圧縮，R特有の型を保持	R専用

3.11.2 分析結果の保存

modelsummary()のoutput引数でファイル形式を指定できる．

# modelsummaryで表をファイルに出力
modelsummary(
  list("Model 1" = model1, "Model 2" = model2),
  output = here("tables", "regression_results.html")
)

# Word形式で出力
modelsummary(
  list("Model 1" = model1, "Model 2" = model2),
  output = here("tables", "regression_results.docx")
)

3.11.3 図の保存

ggsave()で様々な形式で保存できる．プレゼン用はPNG，論文投稿用はPDFが一般的である．

# 図を作成
p <- d |>
  drop_na(f_edu) |>
  ggplot(aes(x = f_edu, y = vocab_w2_g3)) +
  geom_boxplot(fill = "lightblue") +
  labs(x = "父親の学歴", y = "語彙力スコア")

# PNG形式で保存
ggsave(here("figures", "vocab_by_fedu.png"), p,
       width = 6, height = 4, dpi = 300)

# PDF形式で保存（論文向け）
ggsave(here("figures", "vocab_by_fedu.pdf"), p,
       width = 6, height = 4)

ggsave()のポイント

• widthとheightはインチ単位
• dpi = 300で印刷品質（Web用はdpi = 150程度でOK）
• 最後に作成した図はpを省略可能（ggsave("file.png")）

4 まとめ

4.1 本日学んだこと

1. 準備編
   • RStudioの基本操作
   • プロジェクト管理
   • データの読み込み
2. 基礎編
   • dplyrによるデータ操作（filter, select, mutate, summarise）
   • case_match()とcase_when()による変数のリコード
   • 自作関数の作成
   • 記述統計の算出
   • ggplot2による可視化
3. 実践編
   • クロス表分析とカイ二乗検定
   • 回帰分析（単回帰・重回帰・交互作用）
   • ロジスティック回帰と限界効果
   • modelsummaryによるモデル比較
   • パネルデータの作成と分析
   • データの保存・エクスポート

4.2 発展的内容

本コースで扱いきれなかった発展的なトピックを紹介する．

4.2.1 サンプリングウェイトの構成

調査データに付与されるウェイトは、一般に以下の3つの要素の積で構成される。

\[w = w_b \times w_r \times w_c\]

記号	名称	説明
\(w_b\)	デザインウェイト	標本設計により決定。各標本の抽出確率の逆数
\(w_r\)	回答率調整ウェイト	調査拒否・未回答による偏りを調整
\(w_c\)	キャリブレーション・ファクタ	既知の母集団情報（人口総数等）との整合性を調整

デザインウェイトは、二段抽出の場合、層 \(i\) の調査区 \(j\) に属する世帯のウェイトは以下で計算される：

\[w_{b_{ij}} = \frac{M_i}{m_i} \times \frac{N_{ij}}{n_{ij}}\]

ここで、\(M_i\) は層 \(i\) の総調査区数、\(m_i\) は標本調査区数、\(N_{ij}\) は調査区 \(ij\) の総世帯数、\(n_{ij}\) は標本世帯数である。

回答率調整ウェイトは、標本数 \(n_{ij}\) に対して回答数が \(n'_{ij}\) の場合：

\[w_{r_{ij}} = \frac{n_{ij}}{n'_{ij}}\]

キャリブレーション・ファクタは、人口予測値等のベンチマークを用いてウェイトを調整する。

ウェイトの提供状況

公開されるミクロデータには、これらを統合した最終ウェイトのみが提供されることが多い。ウェイトの構成要素が個別に提供されることは稀であり、ウェイトの計算方法は調査の技術資料を参照する必要がある。

4.2.2 ブートストラップ法とジャックナイフ法

標準誤差や信頼区間の推定に、リサンプリング法が有用である。特に、解析的な標準誤差の計算が困難な場合や、複雑なサンプリングデザインを考慮する場合に用いられる。

4.2.2.1 ブートストラップ法

元の標本から復元抽出でリサンプルを繰り返し生成し、統計量の分布を推定する方法である。

library(boot)

# ブートストラップ用の関数を定義
# data: データ、indices: リサンプルのインデックス
boot_mean <- function(data, indices) {
  mean(data[indices], na.rm = TRUE)
}

# ブートストラップの実行（1000回）
set.seed(123)
boot_result <- boot(data = d$vocab_w2_g9,
                    statistic = boot_mean,
                    R = 1000)

# 結果の確認
boot_result

# 信頼区間（パーセンタイル法、BCa法など）
boot.ci(boot_result, type = c("perc", "bca"))

4.2.2.2 ジャックナイフ法

標本サイズを \(n\) とするとき、1つの観測値を除外したリサンプル（サイズ \(n-1\)）を \(n\) 組作成し、統計量の分布を推定する方法である。

# ジャックナイフ法による標準誤差の推定
jackknife_se <- function(x) {
  n <- length(x)
  x <- x[!is.na(x)]  # 欠損値を除外
  n <- length(x)

  # 各観測値を1つずつ除外した平均を計算
  jack_means <- numeric(n)
  for (i in 1:n) {
    jack_means[i] <- mean(x[-i])
  }

  # ジャックナイフ標準誤差
  se <- sqrt((n - 1) / n * sum((jack_means - mean(jack_means))^2))
  return(se)
}

# 使用例
jackknife_se(d$vocab_w2_g9)

# 比較：通常の標準誤差
sd(d$vocab_w2_g9, na.rm = TRUE) / sqrt(sum(!is.na(d$vocab_w2_g9)))

ブートストラップ法とジャックナイフ法の比較

特徴	ブートストラップ法	ジャックナイフ法
リサンプル数	任意（通常1000〜10000）	\(n\)（観測数と同じ）
抽出方法	復元抽出	1つを除外
計算量	多い	比較的少ない
適用範囲	広い（複雑な統計量に対応）	滑らかな統計量向け
バイアス補正	BCa法などで対応可能	可能

4.2.2.3 回帰分析でのブートストラップ

回帰係数の標準誤差をブートストラップで推定する例：

# 回帰係数のブートストラップ
boot_coef <- function(data, indices) {
  d_boot <- data[indices, ]
  model <- lm(vocab_w2_g9 ~ f_edu + income__2015, data = d_boot)
  return(coef(model))
}

# 分析用データ（完全ケースのみ）
d_complete <- d |>
  select(vocab_w2_g9, f_edu, income__2015) |>
  na.omit()

# ブートストラップの実行
set.seed(123)
boot_reg <- boot(data = d_complete,
                 statistic = boot_coef,
                 R = 1000)

# 各係数の信頼区間
boot.ci(boot_reg, type = "perc", index = 2)  # f_edu専門・短大の係数
boot.ci(boot_reg, type = "perc", index = 3)  # f_edu大学以上の係数

リサンプリング法の活用場面

• 標準誤差の解析的計算が困難な統計量（中央値、四分位範囲など）
• 非正規分布からの標本
• 複雑なサンプリングデザイン（クラスタ抽出など）
• 小標本での信頼区間推定

4.2.3 欠損値を考慮した分析

パネルデータでは欠損値（脱落・無回答）が避けられない．リストワイズ削除（完全ケース分析）はサンプルサイズの減少とバイアスの原因になる．多重代入法（Multiple Imputation） は，欠損値を複数回代入して不確実性を考慮した分析を行う方法である．

• Amelia II：時系列・横断面データに特化した多重代入法．EMアルゴリズムとブートストラップを使用し，高速に処理できる．パネルデータに適している．
  • パッケージ：Amelia
  • 特徴：時系列構造を考慮，高速，連続変数向け

library(Amelia)
# パネルデータの多重代入（5つのデータセットを生成）
amelia_out <- amelia(d_wide,
                     m = 5,                    # 代入データセット数
                     ts = "year",              # 時間変数
                     cs = "PanelID",           # 個人ID
                     noms = c("sex", "f_edu")) # 名義変数
# 代入済みデータの取り出し
amelia_out$imputations[[1]]  # 1つ目の代入データ

• mice：Multivariate Imputation by Chained Equations．変数ごとに条件付き分布を指定し，連鎖方程式で代入する．柔軟性が高く，カテゴリ変数にも対応．
  • パッケージ：mice
  • 特徴：変数の型に応じた代入方法，柔軟なモデル指定

library(mice)
# 多重代入の実行
mice_out <- mice(d_analysis,
                 m = 5,           # 代入データセット数
                 method = "pmm", # 予測平均マッチング
                 maxit = 10)      # 反復回数
# 代入データで回帰分析を実行し，結果を統合
fit <- with(mice_out, lm(vocab_w2_g3 ~ sex + f_edu))
pooled <- pool(fit)
summary(pooled)

Rubinのルール（Rubin’s rules）

多重代入法では，複数の代入データセット（通常5〜20個）それぞれで分析を行い，その結果をRubinのルールに基づいて統合（プール）する．

• 点推定値：各代入データセットの推定値の平均
• 分散の推定：代入内分散（within-imputation variance）と代入間分散（between-imputation variance）を組み合わせる

miceパッケージのpool()関数はRubinのルールを自動的に適用する．

Ameliaとmiceの使い分け

• Amelia：パネルデータ，連続変数が多い場合，高速処理が必要な場合
• mice：横断データ，カテゴリ変数が多い場合，柔軟なモデル指定が必要な場合

4.2.4 脱落パターンの分析

パネル調査では，調査回を重ねるごとに回答者が減少する（脱落，attrition）．脱落がランダムに生じるのか，それとも特定の属性を持つ人が脱落しやすいのかを確認することは，分析結果のバイアスを評価する上で重要である．

本データには各Waveの回答フラグ（w1回答フラグ～w7回答フラグ）が含まれており，親子それぞれの回答状況を把握できる．

• 1: 親子とも回答（小1〜3含む）
• 2: 親のみ回答
• 3: 子どものみ回答
• 4: 回答なし
• 5: 分析対象外
• 6: 調査対象外（まだ調査に入っていないコーホート等）

# 回答フラグを使って協力パターンを分析
# 元データから回答フラグを読み込む
d_raw <- read_dta(here("data", "raw", "1571", "1571.dta"))

# 回答フラグ変数を選択
d_attrition <- d_raw |>
  select(PanelID, starts_with("w") & ends_with("回答フラグ")) |>
  rename(
    resp_w1 = `w1回答フラグ`,
    resp_w2 = `w2回答フラグ`,
    resp_w3 = `w3回答フラグ`,
    resp_w4 = `w4回答フラグ`,
    resp_w5 = `w5回答フラグ`,
    resp_w6 = `w6回答フラグ`,
    resp_w7 = `w7回答フラグ`
  )

# 各Waveの回答状況の分布（全カテゴリ）
d_attrition |>
  pivot_longer(starts_with("resp_"),
               names_to = "wave",
               values_to = "response") |>
  mutate(
    wave = parse_number(wave),
    response_label = case_match(
      response,
      1 ~ "親子とも回答",
      2 ~ "親のみ回答",
      3 ~ "子のみ回答",
      4 ~ "回答なし",
      5 ~ "分析対象外",
      6 ~ "調査対象外"
    ) |> factor(levels = c("親子とも回答", "親のみ回答",
                           "子のみ回答", "回答なし",
                           "分析対象外", "調査対象外"))
  ) |>
  count(wave, response_label) |>
  ggplot(aes(x = wave, y = n, fill = response_label)) +
  geom_col(position = "stack") +
  scale_x_continuous(breaks = 1:7, labels = paste0("W", 1:7)) +
  scale_fill_viridis_d(option = "C", end = 0.8) +
  labs(x = "Wave", y = "人数", fill = "回答状況",
       title = "各Waveの回答状況の分布（全カテゴリ）")

「調査対象外」を除いた調査対象者のみで回答パターンを確認すると，脱落の推移がより明確になる．

# 調査対象者のみ（調査対象外=6を除く）で割合を表示
d_attrition |>
  pivot_longer(starts_with("resp_"),
               names_to = "wave",
               values_to = "response") |>
  filter(response != 6) |>  # 調査対象外を除く
  mutate(
    wave = parse_number(wave),
    response_label = case_match(
      response,
      1 ~ "親子とも回答",
      2 ~ "親のみ回答",
      3 ~ "子のみ回答",
      4 ~ "回答なし",
      5 ~ "分析対象外"
    ) |> factor(levels = c("親子とも回答", "親のみ回答",
                           "子のみ回答", "回答なし", "分析対象外"))
  ) |>
  count(wave, response_label) |>
  mutate(percent = n / sum(n) * 100, .by = wave) |>
  ggplot(aes(x = wave, y = percent, fill = response_label)) +
  geom_col(position = "stack") +
  geom_text(aes(label = sprintf("%.1f%%", percent)),
            position = position_stack(vjust = 0.5),
            size = 3, color = "white") +
  scale_x_continuous(breaks = 1:7, labels = paste0("W", 1:7)) +
  scale_fill_viridis_d(option = "C", end = 0.8) +
  labs(x = "Wave", y = "割合 (%)", fill = "回答状況",
       title = "各Waveの回答状況の割合")

Wave1で調査対象だった人（回答フラグが1〜5）に限定して，脱落パターンを分析する．「親子とも回答」のケースに限定した分析を行う場合，どの程度サンプルが減少するかを確認する．

# Wave1で調査対象だった人のみを抽出
d_attrition_w1 <- d_attrition |>
  filter(resp_w1 != 6)  # Wave1で調査対象外だった人を除く

# 親子とも回答した回数を計算
d_attrition_w1 <- d_attrition_w1 |>
  mutate(
    # 各Waveで親子とも回答したか（1 = 親子とも回答, 0 = それ以外）
    # 調査対象外(6)の場合も0とする
    both_w1 = if_else(resp_w1 == 1, 1, 0),
    both_w2 = if_else(resp_w2 == 1, 1, 0),
    both_w3 = if_else(resp_w3 == 1, 1, 0),
    both_w4 = if_else(resp_w4 == 1, 1, 0),
    both_w5 = if_else(resp_w5 == 1, 1, 0),
    both_w6 = if_else(resp_w6 == 1, 1, 0),
    both_w7 = if_else(resp_w7 == 1, 1, 0),
    # 親子とも回答した総回数
    n_both = both_w1 + both_w2 + both_w3 + both_w4 +
             both_w5 + both_w6 + both_w7,
    # 調査対象だった回数（調査対象外=6以外の回数）
    n_target = (resp_w1 != 6) + (resp_w2 != 6) + (resp_w3 != 6) +
               (resp_w4 != 6) + (resp_w5 != 6) + (resp_w6 != 6) + (resp_w7 != 6)
  )

# 親子とも回答した回数の分布（Wave1調査対象者のみ）
d_attrition_w1 |>
  count(n_both) |>
  mutate(percent = n / sum(n) * 100)

# A tibble: 8 × 3
  n_both     n percent
   <dbl> <int>   <dbl>
1      0   135   0.745
2      1  2979  16.4  
3      2  2574  14.2  
4      3  2201  12.1  
5      4  2016  11.1  
6      5  1933  10.7  
7      6  2435  13.4  
8      7  3849  21.2  

初回調査（Wave1）の属性によって，脱落率に差があるかを検討する．

# Wave1の属性を結合
d_attrition_w1 <- d_attrition_w1 |>
  left_join(d |> select(PanelID, sex, f_edu), by = "PanelID")

# 父親学歴別の平均回答回数
d_attrition_w1 |>
  drop_na(f_edu) |>
  summarise(
    mean_both = mean(n_both, na.rm = TRUE),
    se = sd(n_both, na.rm = TRUE) / sqrt(n()),
    ll = mean_both - 1.96 * se,
    ul = mean_both + 1.96 * se,
    n = n(),
    .by = f_edu
  ) |>
  ggplot(aes(x = f_edu, y = mean_both, ymin = ll, ymax = ul)) +
  geom_col(fill = "steelblue", alpha = 0.7) +
  geom_errorbar(width = 0.2) +
  labs(x = "父親の学歴", y = "親子とも回答した平均回数（7回中）",
       title = "父親学歴別の平均回答回数（95%信頼区間付き）")

# 脱落と初回特性の関連をロジスティック回帰で検討
# 「全7回親子とも回答」 vs 「それ以外」
d_attrition_w1 <- d_attrition_w1 |>
  mutate(
    complete = if_else(n_both == 7, 1, 0),
    f_edu_factor = factor(f_edu)
  )

# 脱落の予測モデル
attrition_model <- glm(complete ~ sex + f_edu_factor,
                       data = d_attrition_w1,
                       family = binomial)
summary(attrition_model)

Call:
glm(formula = complete ~ sex + f_edu_factor, family = binomial, 
    data = d_attrition_w1)

Coefficients:
                       Estimate Std. Error z value Pr(>|z|)    
(Intercept)            -1.29492    0.04089 -31.671  < 2e-16 ***
sex男子                -0.11382    0.03867  -2.943  0.00325 ** 
f_edu_factor専門・短大  0.24341    0.05984   4.068 4.75e-05 ***
f_edu_factor大学・院    0.24747    0.04473   5.532 3.16e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 16276  on 14908  degrees of freedom
Residual deviance: 16233  on 14905  degrees of freedom
  (3213 個の観測値が欠損のため削除されました)
AIC: 16241

Number of Fisher Scoring iterations: 4

# 結果の表示
tidy(attrition_model, conf.int = TRUE) |>
  filter(term != "(Intercept)") |>
  mutate(
    odds_ratio = exp(estimate),
    or_ll = exp(conf.low),
    or_ul = exp(conf.high)
  ) |>
  select(term, odds_ratio, or_ll, or_ul, p.value)

# A tibble: 3 × 5
  term                   odds_ratio or_ll or_ul      p.value
  <chr>                       <dbl> <dbl> <dbl>        <dbl>
1 sex男子                     0.892 0.827 0.963 0.00325     
2 f_edu_factor専門・短大      1.28  1.13  1.43  0.0000475   
3 f_edu_factor大学・院        1.28  1.17  1.40  0.0000000316

脱落分析の解釈

• オッズ比が1より大きい属性は，完了（全回答）しやすい傾向
• オッズ比が1より小さい属性は，脱落しやすい傾向
• 脱落が特定の属性に偏っている場合，リストワイズ削除による分析結果にバイアスが生じる可能性がある
• 対処法：多重代入法（上記参照），逆確率重み付け（IPW），感度分析

4.2.5 多項ロジスティック回帰

従属変数が3カテゴリ以上の名義変数の場合，多項ロジスティック回帰を使用する．nnetパッケージのmultinom()関数が使いやすい．

library(nnet)

# 希望進学段階を従属変数とした多項ロジスティック回帰
# 参照カテゴリは最初のカテゴリになる
model_multi <- multinom(edu_hope ~ sex + f_edu + income__2015,
                        data = d)

# 結果の確認
summary(model_multi)

# broomで整理
tidy(model_multi, conf.int = TRUE) |>
  mutate(
    RRR = exp(estimate),  # 相対リスク比
    RRR_ll = exp(conf.low),
    RRR_ul = exp(conf.high)
  ) |>
  select(y.level, term, RRR, RRR_ll, RRR_ul, p.value)

多項ロジスティック回帰の解釈

多項ロジスティック回帰の係数は相対リスク比（Relative Risk Ratio: RRR）として解釈される．これは「参照カテゴリに対する各カテゴリの相対的なオッズ」を表す．しかし，RRRの直感的な解釈は難しい場合が多い．

実践的には，限界効果（marginal effects）を計算して「説明変数が1単位変化したときに各カテゴリを選択する確率がどれだけ変化するか」を示す方が解釈しやすい．

# 限界効果の計算
library(marginaleffects)
avg_slopes(model_multi)

4.2.6 潜在クラス分析

観測されない集団の異質性（サブグループ）を分類する手法である．poLCAパッケージを使用する．

library(poLCA)

# 潜在クラス分析用のデータ準備
# カテゴリ変数は1から始まる整数にコードする必要がある
d_lca <- d |>
  select(PanelID,
         item1 = ds0000000244_w1c,  # 授業が楽しい
         item2 = ds0000000245_w1c,  # 友だちとすごすのが楽しい
         item3 = ds0000000246_w1c,  # 尊敬できる先生がいる
         item4 = ds0000000249_w1c,  # 自分のクラスが好き
         item5 = ds0000000250_w1c   # 自分の学校が好き
         ) |>
  mutate(across(item1:item5, ~ recode_missing(.x))) |>
  drop_na()

# モデル式
f <- cbind(item1, item2, item3, item4, item5) ~ 1

# 2クラスモデル
lca2 <- poLCA(f, d_lca, nclass = 2, verbose = FALSE)

# 3クラスモデル
lca3 <- poLCA(f, d_lca, nclass = 3, verbose = FALSE)

# BICでモデル比較（値が小さいほど良い）
c(lca2$bic, lca3$bic)

# クラス所属確率
head(lca2$posterior)

潜在クラス分析の活用

• 異質な集団の発見：回答パターンから異なるタイプの集団を識別
• クラス数の決定：BICやAICを比較して最適なクラス数を選択
• プロファイル分析：各クラスの特徴を記述し，意味のあるラベルを付与

4.2.7 潜在変数モデル

観測されない潜在変数を仮定してモデル化する手法群である．

• 確認的因子分析（CFA）：複数の観測変数から潜在変数（因子）を測定する．尺度の妥当性検証などに使用
• 構造方程式モデル（SEM）：潜在変数間の関係をモデル化する．パス解析と因子分析を統合した枠組み
• 成長曲線モデル：個人の変化の軌跡（切片・傾き）を潜在変数として推定し，その個人差を説明する

パッケージ：lavaan（SEM, CFA, 成長曲線モデル），lme4, nlme（成長曲線モデル）

潜在変数モデルと因果推論

これらの手法は変数間の関連構造を記述するものであり，それ自体が因果効果を推定するわけではない．特に交差遅延モデル（cross-lagged panel model）は，時間不変の交絡を統制できないため，因果推論には適さないことが指摘されている．因果効果の推定には，固定効果モデルや傾向スコア法など，因果推論に適した手法を選択する必要がある．

4.2.8 マルチレベルモデル（階層線形モデル）

パネルデータや学校・地域などにネストされたデータでは，マルチレベルモデル（階層線形モデル，混合効果モデル）が有用である．個人レベルと集団レベルの変動を分離し，集団間の異質性を考慮した推定ができる．

library(lme4)

# ランダム切片モデル：個人ごとに切片が異なる
model_ri <- lmer(vocab ~ grade + sex + (1 | PanelID), data = d_long)

# ランダム傾きモデル：個人ごとに切片と傾きが異なる
model_rs <- lmer(vocab ~ grade + sex + (1 + grade | PanelID), data = d_long)

summary(model_rs)

マルチレベルモデルと固定効果モデルの違い

• マルチレベルモデル：集団効果（ランダム効果）を正規分布から抽出されたものとしてモデル化．集団間の変動に関心がある場合や，集団レベルの説明変数の効果を推定したい場合に適する
• 固定効果モデル：観察されない異質性を完全に統制．集団レベルの効果は推定できないが，内生性への対処として強力

4.2.9 ベイズ推定

ベイズ推定は，事前分布と尤度から事後分布を求める統計的推論の枠組みである．回帰分析，固定効果モデル，マルチレベルモデル，潜在変数モデルなど，様々な分析手法に適用できる．

4.2.9.1 brmsによるベイズ推定

brmsパッケージを使うと，Rの回帰モデル構文でベイズ推定が可能である．

library(brms)

# 通常の回帰分析をベイズで推定
model_bayes <- brm(vocab_w2_g3 ~ sex + f_edu + income__2015,
                   data = d,
                   family = gaussian(),
                   seed = 1234)

# 結果の確認
summary(model_bayes)

# 事後分布の可視化
plot(model_bayes)

# マルチレベルモデルのベイズ推定
model_mlm_bayes <- brm(vocab ~ grade + sex + (1 + grade | PanelID),
                        data = d_long,
                        family = gaussian(),
                        seed = 1234)

4.2.9.2 Stanによる直接記述

より柔軟なモデル構築のためには，Stanコードを直接記述することを推奨する．Stanで書くことで，モデルの構造を明示的に理解でき，brmsでは対応しにくい複雑なモデルも実装できる．

library(rstan)
library(cmdstanr)

# Stanコードの例：線形回帰
stan_code <- "
data {
  int<lower=0> N;           // サンプルサイズ
  int<lower=0> K;           // 説明変数の数
  matrix[N, K] X;           // 説明変数行列
  vector[N] y;              // 被説明変数
}

parameters {
  vector[K] beta;           // 回帰係数
  real<lower=0> sigma;      // 残差標準偏差
}

model {
  // 事前分布
  beta ~ normal(0, 10);
  sigma ~ cauchy(0, 5);

  // 尤度
  y ~ normal(X * beta, sigma);
}
"

# コンパイルと推定
# stan_model <- stan_model(model_code = stan_code)
# fit <- sampling(stan_model, data = stan_data)

brmsとStanの使い分け

• brms：Rの回帰モデル構文でStanモデルを記述できる．初学者にも扱いやすく，多くの分析に対応．内部でStanコードを生成するため，stancode(model)で生成されたコードを確認できる
• Stan（rstan/cmdstanr）：Stanコードを直接記述する．モデルの構造を完全に制御でき，より複雑なモデル（カスタム尤度，複雑な階層構造，潜在変数モデルなど）に対応可能

brmsで始めて，生成されたStanコードを読むことでStanの文法を学び，必要に応じてStanで直接書くようにするとよい．長期的にはStanで書けるようになることを推奨する．モデルの構造を深く理解でき，研究の幅が広がる．

ベイズ推定の利点

• 事前知識の活用：先行研究の知見を事前分布として組み込める
• 不確実性の表現：点推定だけでなく，パラメータの事後分布全体を得られる
• 小サンプルでの安定性：事前分布による正則化効果で，小サンプルでも安定した推定が可能
• 複雑なモデルへの対応：最尤法では推定困難なモデルも扱える
• マルチレベルモデルとの相性：グループ数が少ない場合や分散成分の推定が不安定な場合にも頑健

4.2.10 因果推論の基礎

因果推論の枠組みでは，潜在アウトカム（potential outcomes）を用いて因果効果を定義する．個人\(i\)が処置を受けた場合のアウトカム\(Y_i(1)\)と受けなかった場合のアウトカム\(Y_i(0)\)の差\(Y_i(1) - Y_i(0)\)が個人レベルの因果効果である．しかし，同一個人について両方のアウトカムを同時に観測することはできない（因果推論の根本問題）．

識別と推定の区別

因果推論では識別（identification）と推定（estimation）を区別することが重要である．

• 識別：観測データから因果効果を定義できるか．どのような仮定のもとで，潜在アウトカムの比較が観測データの比較に置き換えられるか
• 推定：識別された量をデータからどのように計算するか

回帰分析を実行すれば何らかの係数は得られるが，それが因果効果として解釈できるかは別問題である．推定すれば因果分析ができるわけではない．因果効果の識別には，交絡がないこと（条件付き独立性），処置の割り当てが確率的であること（正値性）などの仮定が必要であり，これらの仮定の妥当性は研究者が実質的な知識に基づいて判断しなければならない．

主な識別戦略

観測データから因果効果を識別するための代表的なアプローチ：

• 選択オブザーバブル：観測された共変量で条件付けることで交絡を除去できると仮定（回帰調整，マッチング，IPW）
• 固定効果：時間不変の交絡因子を個人固定効果として除去
• 操作変数：処置には影響するがアウトカムには直接影響しない変数を利用
• 回帰不連続デザイン：処置割り当ての閾値付近での比較
• 差の差法：処置群と統制群の時間変化の差を比較

いずれの戦略も，特定の仮定のもとでのみ因果効果を識別できる．仮定の妥当性はデータからは検証できず，研究者の判断に委ねられる．

• 回帰代入モデル：アウトカムの条件付き期待値をモデル化し，反実仮想を予測する方法
  • パッケージ：marginaleffects
• 逆確率重み付け（IPW）：処置を受ける確率（傾向スコア）の逆数で重み付けし，交絡を調整する方法
  • パッケージ：WeightIt
• 二重にロバストな方法（AIPW）：回帰モデルとIPWを組み合わせ，どちらか一方が正しければ一致推定量となる方法
  • パッケージ：AIPW, DoubleML

4.2.11 因果推論の発展

• 媒介分析：処置が媒介変数を通じてアウトカムに影響する経路（間接効果）を分解する分析．
  • パッケージ：mediation, CMAverse, rwrmed
• 複数の処置がある変数：3カテゴリ以上の処置変数や，連続的な処置変数を扱う因果推論の手法．
  • パッケージ：ltmle

4.2.12 機械学習による予測

• SuperLearner：複数の機械学習アルゴリズムを組み合わせて最適な予測モデルを構築するアンサンブル学習法．因果推論との組み合わせ（TMLE等）にも活用される．
  • パッケージ：SuperLearner, sl3

4.2.13 研究の透明性と再現性

二次分析においても，研究の透明性と再現可能性を高める取り組みが求められている．

• データ管理計画（DMP: Data Management Plan）：研究データの収集，管理，共有，保存に関する計画を事前に文書化する．研究開始前にデータの取り扱い方針を明確にすることで，分析の一貫性と再現性を担保する．
  • 二次分析では，使用するデータの出典，変数の選択基準，前処理の手順を明記することが重要
  • 日本学術振興会の科研費申請でもDMPの提出が求められる場合がある
• 事前登録（Preregistration）：分析計画を事前に登録し，研究の仮説・分析手法を固定する．結果を見てから仮説を後付けするHARKing（Hypothesizing After the Results are Known）を防ぐ．
  • プラットフォーム：OSF，AsPredicted
  • 二次分析でも，どの変数をどのように分析するかを事前に登録することで，探索的分析と確証的分析を区別できる

事前登録の実践

二次分析で事前登録を行う場合，以下の項目を含めるとよい：

1. 使用するデータセットと対象年次
2. 分析対象のサブサンプル（学年，性別等の限定条件）
3. 主要な変数の定義（どの変数をどう加工するか）
4. 検定する仮説と統計モデル
5. 有意水準と多重比較の調整方法

4.2.14 パッケージの中身から学ぶ

Rの関数がどのように実装されているかを確認することで，統計手法の理解を深めることができる．

# 関数名を()なしで入力すると中身が表示される
mean

UseMethod("mean")のように表示される場合は，methods()で具体的なメソッドを確認する．

# 利用可能なメソッドを確認
methods(mean)

# デフォルトのメソッドを確認
mean.default

# getAnywhere()でも確認できる
getAnywhere(mean.default)

パッケージ内の非公開関数（エクスポートされていない関数）は:::で参照できる．

# パッケージ名:::関数名 で非公開関数を参照
# 例：DescToolsパッケージのOddsRatio関数の中身
DescTools:::OddsRatio.default

関数の中身を読むメリット

• 統計手法の計算過程を具体的に理解できる
• デフォルトの設定値や前提条件を確認できる
• 自分で関数を作成する際の参考になる
• パッケージのバグや想定外の挙動を発見できることがある

4.3 参考資料

4.3.1 書籍

• 松村優哉ほか（2021）『改訂2版 RユーザのためのRStudio[実践]入門』技術評論社
• Hadley Wickham・Garrett Grolemund（2017）『Rではじめるデータサイエンス』オライリー・ジャパン

4.3.2 オンライン

• R for Data Science (2e)：tidyverseを使ったデータ分析の入門書（無料）
• Quarto：Quarto公式ドキュメント
• ggplot2公式サイト：関数リファレンスと作例
• CRAN Task Views：分野別のパッケージ一覧

4.4 データの利用について

• 本講義で使用したSSJデータアーカイブのデータは，3月中に削除すること
• 研究目的での継続利用を希望される場合は，SSJデータアーカイブに申請すること
• https://csrda.iss.u-tokyo.ac.jp/

4.5 寄付について

• 社会科学データサイエンスの基盤構築基金
• https://utf.u-tokyo.ac.jp/project/pjt202

5 付録A：各Waveの職業階層変数の作成

本文では2015年（W1）の職業・就業形態のみを用いたが、パネルデータの特性を活かすために各Waveの職業階層変数を作成する方法を示す。

5.1 各Waveの変数作成

本文と同じリコード処理を各年度に適用する。変数名は年度ベース（例：f_occ__2016、f_class__2016）とする。

# 2016年（W2）の父親の職種
d <- d |>
  mutate(
    f_occ__2016 = case_match(
      zap_labels(PFworkjob16),
      1 ~ "事務・営業",
      2 ~ "販売・サービス",
      3 ~ "技能・労務",
      4 ~ "保安",
      5 ~ "運輸",
      6 ~ "農林漁業",
      7 ~ "専門・技術",
      8 ~ "管理職",
      c(9, 10) ~ NA,
      .default = NA
    ),
    f_occ__2016 = factor(f_occ__2016)
  )

# 2016年（W2）の父親の就業形態
d <- d |>
  mutate(
    f_emp__2016 = case_match(
      zap_labels(PFworktype16),
      1 ~ "正規雇用",
      c(2, 3, 4) ~ "非正規雇用",
      5 ~ "自営業",
      c(6, 8) ~ NA,
      7 ~ "無職",
      .default = NA
    ),
    f_emp__2016 = if_else(zap_labels(PFworktype16) == 8888, "父親なし", f_emp__2016),
    f_emp__2016 = factor(f_emp__2016)
  )

# 2016年（W2）の父親の職業階層
d <- d |>
  mutate(
    f_class__2016 = case_when(
      f_occ__2016 == "管理職" & f_emp__2016 %in% c("正規雇用", "自営業") ~ "専門・管理職",
      f_occ__2016 == "専門・技術" & f_emp__2016 %in% c("正規雇用", "自営業") ~ "専門・管理職",
      f_occ__2016 == "事務・営業" & f_emp__2016 == "正規雇用" ~ "事務・販売職",
      f_occ__2016 == "販売・サービス" & f_emp__2016 == "正規雇用" ~ "事務・販売職",
      f_emp__2016 == "自営業" ~ "自営業",
      f_occ__2016 == "農林漁業" & f_emp__2016 == "自営業" ~ "自営業",
      f_occ__2016 %in% c("技能・労務", "運輸", "保安") & f_emp__2016 == "正規雇用" ~ "マニュアル",
      f_occ__2016 == "農林漁業" & f_emp__2016 == "正規雇用" ~ "マニュアル",
      f_emp__2016 == "非正規雇用" ~ "非正規雇用",
      f_emp__2016 == "無職" ~ "無職",
      f_emp__2016 == "父親なし" ~ "父親なし",
      .default = NA
    ),
    f_class__2016 = factor(f_class__2016, levels = c(
      "専門・管理職", "事務・販売職", "自営業",
      "マニュアル", "非正規雇用", "無職", "父親なし"
    ))
  )

d |>
  count(f_class__2016) |>
  mutate(pct = n / sum(n) * 100)

# A tibble: 8 × 3
  f_class__2016     n    pct
  <fct>         <int>  <dbl>
1 専門・管理職   4592 15.4  
2 事務・販売職   4704 15.8  
3 自営業         1122  3.76 
4 マニュアル     3037 10.2  
5 非正規雇用      366  1.23 
6 無職            106  0.355
7 父親なし        945  3.17 
8 <NA>          14974 50.2  

同様に、他の年度（2017〜2021年）や母親についても作成できる。

5.2 職業階層の変化の検出

本文で作成したf_class（2015年）と付録で作成したf_class__2016を使って、職業階層の変化を検出できる。

# 2015年から2016年への変化
d <- d |>
  mutate(
    f_class_changed = f_class != f_class__2016 &
                      !is.na(f_class) & !is.na(f_class__2016)
  )

# 変化した件数と割合
d |>
  filter(!is.na(f_class) & !is.na(f_class__2016)) |>
  count(f_class_changed) |>
  mutate(pct = n / sum(n) * 100)

# A tibble: 2 × 3
  f_class_changed     n   pct
  <lgl>           <int> <dbl>
1 FALSE            8053  76.9
2 TRUE             2424  23.1

# どの階層間の移動が多いか（件数と割合）
d |>
  filter(f_class_changed == TRUE) |>
  count(f_class, f_class__2016, sort = TRUE) |>
  mutate(pct = n / sum(n) * 100) |>
  head(10)

# A tibble: 10 × 4
   f_class      f_class__2016     n   pct
   <fct>        <fct>         <int> <dbl>
 1 事務・販売職 専門・管理職    383 15.8 
 2 専門・管理職 事務・販売職    371 15.3 
 3 専門・管理職 マニュアル      331 13.7 
 4 マニュアル   専門・管理職    309 12.7 
 5 専門・管理職 自営業          132  5.45
 6 自営業       専門・管理職    131  5.40
 7 事務・販売職 マニュアル      124  5.12
 8 マニュアル   事務・販売職    113  4.66
 9 専門・管理職 父親なし         58  2.39
10 事務・販売職 父親なし         55  2.27

職業階層の変化を用いた分析例

• 親の失業（マニュアル→無職）が子どもの学力に与える影響
• 母親の就業開始（無職→非正規雇用）と子どもの生活時間の変化
• 親の昇進（マニュアル→専門・管理職）と教育投資の変化

6 付録B：中学受験に関する変数

パネルデータの特性を活かした分析例として、中学受験に関する変数を取り上げる。

6.1 中学受験関連変数の確認

このデータセットには、保護者調査で中学受験に関する質問が含まれている。

• ds0000000187_w*p：中学受験予定（小学生の保護者への質問）
• ds0000000181_w*p：中学受験有無（中学生の保護者への質問）

分析対象の限定

中学受験に関する質問は、対象となる学年が限られている。

• 中学受験予定：小学1〜6年生（学年1-6）の保護者に質問
• 中学受験有無：中学1〜3年生（学年7-9）の保護者に質問

したがって、分析対象は限定されたサンプルとなる点に注意が必要である。

# W1（2015年）の中学受験有無（中学生の保護者）
# as_factor()で値ラベルをそのままfactorに変換し、不要なカテゴリをNAに
d <- d |>
  mutate(
    juken_result__2015 = as_factor(ds0000000181_w1p),
    juken_result__2015 = na_if(juken_result__2015, "質問なし"),
    juken_result__2015 = na_if(juken_result__2015, "非該当"),
    juken_result__2015 = na_if(juken_result__2015, "無回答・不明"),
    juken_result__2015 = fct_drop(juken_result__2015)
  )

# 学年7-9（中学生）に限定して確認
d |>
  filter(`w1学年` %in% 7:9) |>
  count(juken_result__2015) |>
  mutate(pct = n / sum(n) * 100)

# A tibble: 5 × 3
  juken_result__2015           n    pct
  <fct>                    <int>  <dbl>
1 中学受験した               752 16.6  
2 準備をしたが途中でやめた    29  0.642
3 考えたがしなかった         294  6.50 
4 考えなかった              3010 66.6  
5 <NA>                       435  9.62 

6.2 全ウェーブを使った出生コーホート別の中学受験率

各ウェーブの中学受験有無を使って、出生コーホート（出生年度）別に中学受験率を算出する。中学1年生になる年度はbirth_year + 13で計算できる。

中学受験有無は中学生の保護者に質問されるため、W1時点で中学生だった2000-2002年生まれはW1で、それ以降の出生コーホートは中1になった時点のウェーブでデータを取得できる。

出生年度	中1になる年度	該当ウェーブ	備考
2000	2013	W1	W1時点で中3
2001	2014	W1	W1時点で中2
2002	2015	W1	W1時点で中1
2003	2016	W2
2004	2017	W3
2005	2018	W4
2006	2019	W5
2007	2020	W6
2008	2021	W7

# 各ウェーブの中学受験有無を変換
d <- d |>
  mutate(
    # W2-W7の中学受験有無をas_factorで変換
    juken_result__2016 = as_factor(ds0000000181_w2p),
    juken_result__2017 = as_factor(ds0000000181_w3p),
    juken_result__2018 = as_factor(ds0000000181_w4p),
    juken_result__2019 = as_factor(ds0000000181_w5p),
    juken_result__2020 = as_factor(ds0000000181_w6p),
    juken_result__2021 = as_factor(ds0000000181_w7p)
  )

# 出生コーホート別に中学受験有無を統合
# 中1になる年度（birth_year + 13）に応じて該当するウェーブの値を取得
d <- d |>
  mutate(
    # 中1になる年度
    jhs_entry_year = birth_year + 13,
    # 該当ウェーブの中学受験有無を取得
    juken_result = case_when(
      jhs_entry_year <= 2015 ~ juken_result__2015,
      jhs_entry_year == 2016 ~ juken_result__2016,
      jhs_entry_year == 2017 ~ juken_result__2017,
      jhs_entry_year == 2018 ~ juken_result__2018,
      jhs_entry_year == 2019 ~ juken_result__2019,
      jhs_entry_year == 2020 ~ juken_result__2020,
      jhs_entry_year == 2021 ~ juken_result__2021,
      .default = NA
    ),
    # 中学受験したかどうかの二値変数
    juken_did = case_when(
      juken_result == "中学受験した" ~ 1,
      juken_result == "中学受験をした" ~ 1,
      juken_result %in% c("準備をしたが途中でやめた", "中学受験の準備をしたが、途中でやめた",
                          "考えたがしなかった", "中学受験を考えたが、しなかった",
                          "考えなかった", "中学受験は考えなかった") ~ 0,
      .default = NA
    )
  )

# 出生コーホート別の中学受験率
d |>
  drop_na(juken_did, birth_year) |>
  group_by(birth_year) |>
  summarise(
    n = sum(!is.na(juken_did)),
    juken_rate = mean(juken_did, na.rm = TRUE) * 100
  )

# A tibble: 10 × 3
   birth_year     n juken_rate
        <dbl> <int>      <dbl>
 1       1999     1        0  
 2       2000  1379       17.3
 3       2001  1369       19.4
 4       2002  1337       18.5
 5       2003  1125       18.3
 6       2004  1069       15.2
 7       2005   978       17.9
 8       2006  1109       16.1
 9       2007  1144       17.7
10       2008  1266       16.8

# 父親の学歴別の中学受験率
d |>
  drop_na(juken_did, f_edu) |>
  group_by(f_edu) |>
  summarise(
    n = sum(!is.na(juken_did)),
    juken_rate = mean(juken_did, na.rm = TRUE) * 100
  )

# A tibble: 3 × 3
  f_edu          n juken_rate
  <fct>      <int>      <dbl>
1 中学・高校  2795       8.59
2 専門・短大  1543      11.3 
3 大学・院    4991      24.7 

# 父親の職業階層別の中学受験率
d |>
  drop_na(juken_did, f_class) |>
  group_by(f_class) |>
  summarise(
    n = sum(!is.na(juken_did)),
    juken_rate = mean(juken_did, na.rm = TRUE) * 100
  )

# A tibble: 6 × 3
  f_class          n juken_rate
  <fct>        <int>      <dbl>
1 専門・管理職  2971      22.8 
2 事務・販売職  3006      18.8 
3 自営業         741      14.6 
4 マニュアル    1893       9.61
5 非正規雇用     223       7.62
6 無職            65      15.4 

6.3 中学受験率の可視化

# 出生コーホート別・父親の学歴別の中学受験率（信頼区間付き）
d |>
  drop_na(juken_did, f_edu, birth_year) |>
  filter(birth_year >= 2000 & birth_year <= 2008) |>
  group_by(birth_year, f_edu) |>
  summarise(
    n = sum(!is.na(juken_did)),
    juken_rate = mean(juken_did, na.rm = TRUE),
    se = sqrt(juken_rate * (1 - juken_rate) / n),
    ci_lower = juken_rate - 1.96 * se,
    ci_upper = juken_rate + 1.96 * se,
    .groups = "drop"
  ) |>
  mutate(
    juken_rate = juken_rate * 100,
    ci_lower = ci_lower * 100,
    ci_upper = ci_upper * 100
  ) |>
  ggplot(aes(x = birth_year, y = juken_rate, color = f_edu,
             shape = f_edu, linetype = f_edu)) +
  geom_line() +
  geom_point(size = 2) +
  geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper, fill = f_edu),
              alpha = 0.2, color = NA) +
  scale_x_continuous(breaks = 2000:2008) +
  scale_color_viridis_d(option = "C", end = 0.8) +
  scale_fill_viridis_d(option = "C", end = 0.8) +
  labs(
    x = "出生年度",
    y = "中学受験率（%）",
    color = "父親の学歴",
    fill = "父親の学歴",
    shape = "父親の学歴",
    linetype = "父親の学歴",
    title = "出生コーホート別・父親の学歴別の中学受験率",
    subtitle = "95%信頼区間付き"
  ) +
  ylim(0, 50)

中学受験を用いた発展的な分析例

• 中学受験経験と中学校での学力・学習時間の関連
• 受験予定の変化（やめた/新たに始めた）と家庭環境の変化
• 中学受験と教育費支出の関連
• 兄弟姉妹の有無と中学受験率

6.4 付録C：家族構成に関する変数

配偶者の有無、子ども人数、出生順位など、家族構成に関する変数を作成する。

• 時間不変の変数：出生順位、第一子かどうか（調査対象児の属性）
• 時変の変数：配偶者の有無、子どもの人数、一人っ子かどうか（家族状況は変化しうる）

6.4.1 出生順位（時間不変）

出生順位は調査対象児の属性であり、時間によって変化しない。W2以降に追加されたコーホートもあるため、coalesce()で最初に得られた値を採用する。

# 各Waveの出生順位を確認（W1）
d |>
  count(ds0000000079_w1p) |>
  mutate(pct = n / sum(n) * 100)

# A tibble: 10 × 3
   ds0000000079_w1p        n      pct
   <dbl+lbl>           <int>    <dbl>
 1    1 [第１子]        8557 28.7    
 2    2 [第２子]        6305 21.1    
 3    3 [第３子]        1635  5.48   
 4    4 [第４子]         207  0.694  
 5    5 [第５子]          17  0.0570 
 6    6 [第６子]           6  0.0201 
 7    7 [第７子]           2  0.00670
 8   10 [第10子]           1  0.00335
 9 9999 [無回答・不明]    80  0.268  
10   NA                13036 43.7    

# 時間不変の出生順位変数を作成
# coalesce()で最初のWaveの値を優先的に採用
d <- d |>
  mutate(
    # 各Waveの出生順位（欠損値処理）
    birth_order_w1 = na_if(ds0000000079_w1p, 7777) |> na_if(9999),
    birth_order_w2 = na_if(ds0000000079_w2p, 7777) |> na_if(9999),
    birth_order_w3 = na_if(ds0000000079_w3p, 7777) |> na_if(9999),
    birth_order_w4 = na_if(ds0000000079_w4p, 7777) |> na_if(9999),
    birth_order_w5 = na_if(ds0000000079_w5p, 7777) |> na_if(9999),
    birth_order_w6 = na_if(ds0000000079_w6p, 7777) |> na_if(9999),
    birth_order_w7 = na_if(ds0000000079_w7p, 7777) |> na_if(9999),
    # 最初のWaveの値を採用
    birth_order = coalesce(birth_order_w1, birth_order_w2, birth_order_w3,
                           birth_order_w4, birth_order_w5, birth_order_w6, birth_order_w7)
  ) |>
  select(-starts_with("birth_order_w"))  # 作業用変数を削除

# 確認
d |>
  count(birth_order) |>
  mutate(pct = n / sum(n) * 100)

# A tibble: 9 × 3
  birth_order     n      pct
  <dbl+lbl>   <int>    <dbl>
1  1 [第１子] 13984 46.9    
2  2 [第２子] 10064 33.7    
3  3 [第３子]  2730  9.15   
4  4 [第４子]   355  1.19   
5  5 [第５子]    49  0.164  
6  6 [第６子]    12  0.0402 
7  7 [第７子]     3  0.0101 
8 10 [第10子]     1  0.00335
9 NA           2648  8.87   

6.4.2 第一子かどうか（時間不変）

# 第一子フラグを作成（出生順位から派生）
d <- d |>
  mutate(
    first_child = case_when(
      birth_order == 1 ~ 1,
      birth_order >= 2 ~ 0,
      .default = NA
    )
  )

# 確認
d |>
  count(first_child) |>
  mutate(pct = n / sum(n) * 100)

# A tibble: 3 × 3
  first_child     n   pct
        <dbl> <int> <dbl>
1           0 13214 44.3 
2           1 13984 46.9 
3          NA  2648  8.87

6.4.3 配偶者の有無（時変）

配偶者の有無は離婚・再婚などにより変化しうるため、各Waveで変数を作成する。

# 配偶者の有無（各Wave）
d <- d |>
  mutate(
    spouse__2015 = case_match(ds0000000025_w1p, 1 ~ 1, 2 ~ 0, .default = NA),
    spouse__2016 = case_match(ds0000000025_w2p, 1 ~ 1, 2 ~ 0, .default = NA),
    spouse__2017 = case_match(ds0000000025_w3p, 1 ~ 1, 2 ~ 0, .default = NA),
    spouse__2018 = case_match(ds0000000025_w4p, 1 ~ 1, 2 ~ 0, .default = NA),
    spouse__2019 = case_match(ds0000000025_w5p, 1 ~ 1, 2 ~ 0, .default = NA),
    spouse__2020 = case_match(ds0000000025_w6p, 1 ~ 1, 2 ~ 0, .default = NA),
    spouse__2021 = case_match(ds0000000025_w7p, 1 ~ 1, 2 ~ 0, .default = NA)
  )

# W1の確認
d |>
  count(spouse__2015) |>
  mutate(pct = n / sum(n) * 100)

# A tibble: 3 × 3
  spouse__2015     n   pct
         <dbl> <int> <dbl>
1            0   893  2.99
2            1 15858 53.1 
3           NA 13095 43.9 

6.4.4 子どもの人数（時変）

子どもの人数は増加しうるため、各Waveで変数を作成する。

# 子どもの人数（各Wave）
d <- d |>
  mutate(
    n_children__2015 = na_if(ds0000000057_w1p, 7777) |> na_if(9999),
    n_children__2016 = na_if(ds0000000057_w2p, 7777) |> na_if(9999),
    n_children__2017 = na_if(ds0000000057_w3p, 7777) |> na_if(9999),
    n_children__2018 = na_if(ds0000000057_w4p, 7777) |> na_if(9999),
    n_children__2019 = na_if(ds0000000057_w5p, 7777) |> na_if(9999),
    n_children__2020 = na_if(ds0000000057_w6p, 7777) |> na_if(9999),
    n_children__2021 = na_if(ds0000000057_w7p, 7777) |> na_if(9999)
  )

# W1の確認
d |>
  count(n_children__2015) |>
  mutate(pct = n / sum(n) * 100)

# A tibble: 10 × 3
   n_children__2015     n      pct
   <dbl+lbl>        <int>    <dbl>
 1  1                2162  7.24   
 2  2                8868 29.7    
 3  3                4588 15.4    
 4  4                 797  2.67   
 5  5                 159  0.533  
 6  6                  41  0.137  
 7  7                  18  0.0603 
 8  8                   3  0.0101 
 9  9                   1  0.00335
10 NA               13209 44.3    

6.4.5 一人っ子かどうか（時変）

一人っ子かどうかは子どもの人数から派生し、弟妹の誕生により変化しうる。

# 一人っ子フラグ（各Wave）
d <- d |>
  mutate(
    only_child__2015 = case_when(n_children__2015 == 1 ~ 1, n_children__2015 >= 2 ~ 0, .default = NA),
    only_child__2016 = case_when(n_children__2016 == 1 ~ 1, n_children__2016 >= 2 ~ 0, .default = NA),
    only_child__2017 = case_when(n_children__2017 == 1 ~ 1, n_children__2017 >= 2 ~ 0, .default = NA),
    only_child__2018 = case_when(n_children__2018 == 1 ~ 1, n_children__2018 >= 2 ~ 0, .default = NA),
    only_child__2019 = case_when(n_children__2019 == 1 ~ 1, n_children__2019 >= 2 ~ 0, .default = NA),
    only_child__2020 = case_when(n_children__2020 == 1 ~ 1, n_children__2020 >= 2 ~ 0, .default = NA),
    only_child__2021 = case_when(n_children__2021 == 1 ~ 1, n_children__2021 >= 2 ~ 0, .default = NA)
  )

# W1の確認
d |>
  count(only_child__2015) |>
  mutate(pct = n / sum(n) * 100)

# A tibble: 3 × 3
  only_child__2015     n   pct
             <dbl> <int> <dbl>
1                0 14475 48.5 
2                1  2162  7.24
3               NA 13209 44.3 

家族構成変数を用いた分析例

• 一人っ子と複数きょうだいでの教育投資の違い
• 出生順位と学力・進学希望の関連
• きょうだい数と通塾率の関連
• ひとり親世帯と両親世帯での教育格差
• 弟妹の誕生が学業成績に与える影響（固定効果モデル）

7 付録D：メタデータの管理とコードブック作成

Rでは、SPSSやStataと異なり、変数ラベルや値ラベルがデータに埋め込まれていない。そのため、メタデータの管理とコードブック作成は分析者自身が行う必要がある。

7.1 変数ラベルと値ラベルの付与

labelledパッケージを使用して、変数にラベルを付与できる。

library(labelled)

# 変数ラベルの付与
var_label(d$vocab_w1_g7) <- "語彙スコア（中1）"
var_label(d$juku__2015) <- "通塾の有無（2015年）"

# 複数変数に一括で付与
var_label(d) <- list(
  vocab_w1_g7 = "語彙スコア（中1）",
  vocab_w2_g9 = "語彙スコア（中3）",
  juku__2015 = "通塾の有無（2015年）",
  f_edu = "父親学歴"
)

# 値ラベルの付与
val_labels(d$juku__2015) <- c("通塾なし" = 0, "通塾あり" = 1)

# 値ラベルの確認
val_labels(d$juku__2015)

havenで読み込んだStataファイルのラベル

haven::read_dta()で読み込んだデータには、すでにStataの変数ラベル・値ラベルがlabelled形式で付与されている。var_label()やval_labels()で確認できる。

7.2 コードブックの生成

7.2.1 sjPlotパッケージ

sjPlot::view_df()は、データフレームの変数一覧をHTML形式で出力する。

library(sjPlot)

# HTMLでコードブックを表示（ブラウザで開く）
view_df(d)

# HTMLファイルとして保存
view_df(d, file = "codebook.html")

7.2.2 codebookパッケージ

より詳細なコードブックを生成する。

library(codebook)

# Rmarkdownベースのコードブック生成
codebook(d)

7.2.3 自作のコードブック関数

シンプルなコードブックを自作することもできる。

# 変数一覧をデータフレームとして作成
make_codebook <- function(data) {
  tibble::tibble(
    変数名 = names(data),
    ラベル = purrr::map_chr(data, ~{
      lbl <- labelled::var_label(.)
      if (is.null(lbl)) NA_character_ else lbl
    }),
    型 = purrr::map_chr(data, ~class(.)[1]),
    欠損数 = purrr::map_int(data, ~sum(is.na(.))),
    ユニーク数 = purrr::map_int(data, ~dplyr::n_distinct(., na.rm = TRUE))
  )
}

# コードブック作成
cb <- make_codebook(d)

# HTMLテーブルとして表示
knitr::kable(cb)

# CSVとして保存
readr::write_csv(cb, "codebook.csv")

7.2.4 skimrパッケージ

変数の記述統計を含む要約を作成する。

library(skimr)

# データの要約
skim(d)

# 特定の変数のみ
d |>
  select(vocab_w1_g7, vocab_w2_g9, juku__2015) |>
  skim()

7.3 値ラベルの一覧表示

# 特定変数の値ラベルを確認
labelled::val_labels(d$f_edu)

# すべてのlabelled変数の値ラベルを取得
labelled::look_for(d)

# キーワードで変数を検索
labelled::look_for(d, "学歴")

look_for()関数

labelled::look_for()はStataのlookforコマンドに相当し、変数名やラベルからキーワード検索ができる。大規模データで特定の変数を探す際に便利。

7.4 メタデータ付きデータの保存

# Stata形式で保存（ラベル情報を保持）
haven::write_dta(d, "data_with_labels.dta")

# SPSS形式で保存
haven::write_sav(d, "data_with_labels.sav")

# RDS形式で保存（R専用だがすべての属性を保持）
saveRDS(d, "data_with_labels.rds")

CSV保存時の注意

write_csv()でCSV形式に保存すると、変数ラベル・値ラベルの情報は失われる。メタデータを保持したい場合は、Stata/SPSS形式かRDS形式で保存する。

7.5 データクリーニングの記録

分析の再現性のため、データクリーニングの手順はスクリプトとして保存する。

# cleaning.R の例
library(tidyverse)
library(haven)
library(labelled)

# 元データ読み込み
d_raw <- read_dta("data/raw/original.dta")

# クリーニング処理
d <- d_raw |>
  # 欠損値コードをNAに変換
  mutate(across(where(is.numeric), ~na_if(., 9999))) |>
  mutate(across(where(is.numeric), ~na_if(., 8888))) |>
  mutate(across(where(is.numeric), ~na_if(., 7777))) |>
  # 新変数の作成

  mutate(
    f_edu = case_match(
      father_education,
      1:2 ~ "中学・高校",
      3:4 ~ "専門・短大",
      5:6 ~ "大学以上"
    ) |> factor(levels = c("中学・高校", "専門・短大", "大学以上"))
  )

# 変数ラベルの付与
var_label(d$f_edu) <- "父親学歴（3区分）"

# クリーニング済みデータの保存
write_dta(d, "data/processed/cleaned.dta")
saveRDS(d, "data/processed/cleaned.rds")

データクリーニングのベストプラクティス

1. 元データは絶対に上書きしない（data/raw/に保管）
2. クリーニング手順はすべてスクリプトに記録
3. クリーニング済みデータは別フォルダに保存（data/processed/）
4. 変数の作成・変換の理由をコメントで記録
5. コードブックも一緒に更新する